## [1] "2019-01-27-17-08-53"

Working Directory

## [1] "/home/peterslab/Alex Alleman/Statewide Microbiome Analysis/Statewide analysis"
set.seed(8765)

Load packages

library(ggplot2)
library(data.table)
library(vegan)
library(ggvegan)
library(dplyr)
library(scales)
library(grid)
library(reshape2)
library(phyloseq)
library(ggpubr)
library(RColorBrewer)
library(ape)
library(grid)
library(knitr)
library(igraph)
library(Matrix)
library(ggnetwork)
library(intergraph)
library(parallel)
library(tinytex)

Colors

farm_col<-(c("#8c510a", "#d8b365", "#f6e8c3", "#f5f5f5", "#c7eae5", "#5ab4ac", "#01665e"))
farm_col_dark<-brewer.pal(7, "Dark2")
farm_col_paired<-brewer.pal(10, "Paired")

Load OTU, Taxa, and Meta data

Add OTU table with sample names on top and OTU names as row names

OTU_nifH<- read.delim(
  "~/Alex Alleman/Statewide Microbiome Analysis/Statewide analysis/nifH_OTUallspring.csv",
  row.names = 1)
head(OTU_nifH)[,1:10]
##      JZ017 JZ018 JZ019 JZ020 JZ021 JZ022 JZ023 JZ024 JZ025 JZ026
## OTU1    21 32923 48199    22   117   103   147  4054  5597   246
## OTU2 35281     5    15  1075   101  1335 26249   376  1705   256
## OTU3    11    14    13    49    20    48     9     8    16     9
## OTU4    80     4    14  2576   144   136   150   129 29299 77807
## OTU5     6     2     7   146   749   526   432  1115 13256     0
## OTU6 19317 11274   799   748  1532   348    32    21    28  2517

These taxa were also created by Mr. DNA through a blast program and the NCBI database

tax_nifH<- read.delim(
  "~/Alex Alleman/Statewide Microbiome Analysis/Statewide analysis/nifH_OTU_ids_2016_fixed1.txt", 
  row.names = 1)
head(tax_nifH)[,1:8]
##          kingdom            phylum                  class          order
## OTU1 k__bacteria     p__firmicutes             c__bacilli  o__bacillales
## OTU2 k__bacteria     p__firmicutes             c__bacilli  o__bacillales
## OTU3 k__bacteria p__proteobacteria c__alphaproteobacteria o__rhizobiales
## OTU4 k__bacteria     p__firmicutes             c__bacilli  o__bacillales
## OTU5 k__bacteria p__proteobacteria c__alphaproteobacteria o__rhizobiales
## OTU6 k__bacteria p__proteobacteria c__alphaproteobacteria o__rhizobiales
##                    family               genus
## OTU1  f__paenibacillaceae    g__paenibacillus
## OTU2  f__paenibacillaceae    g__paenibacillus
## OTU3      f__rhizobiaceae        g__rhizobium
## OTU4  f__paenibacillaceae    g__paenibacillus
## OTU5 f__bradyrhizobiaceae   g__bradyrhizobium
## OTU6 f__bradyrhizobiaceae g__rhodopseudomonas
##                            species                       strain
## OTU1     s__paenibacillus graminis       paenibacillus graminis
## OTU2     s__paenibacillus polymyxa       paenibacillus polymyxa
## OTU3    s__rhizobium leguminosarum      rhizobium leguminosarum
## OTU4       s__paenibacillus terrae paenibacillus terrae hpl_003
## OTU5         s__bradyrhizobium sp.   bradyrhizobium sp. stm3074
## OTU6 s__rhodopseudomonas palustris   rhodopseudomonas palustris

Meta data set has be placed together from all the spring and summer data with excel

meta<- read.delim(
  "~/Alex Alleman/Statewide Microbiome Analysis/Statewide analysis/all_metadata_summer.csv",                    row.names = 1, na.strings = "NA", sep = ",",
                  colClasses = c(rep('factor', 7), rep('numeric', 3), rep('factor', 4), 'numeric',
                                 rep('factor', 3), rep('numeric', 28) ) )
head(meta)[,1:5]
##            Site   ARC Season Sample_dates  Pea_variety
## JZ032 Kalispell NWARC Summer  2016-summer        Delta
## JZ031 Kalispell NWARC Summer  2016-summer  CDC Saffron
## JZ030 Kalispell NWARC Summer  2016-summer AC Earlystar
## JZ034 Kalispell NWARC Summer  2016-summer      Majoret
## JZ033 Kalispell NWARC Summer  2016-summer   DS Admiral
## JZ035 Kalispell NWARC Summer  2016-summer      Navarro

Removed all Havre for analysis because samples were taken in 6" sections instead of 12"

meta2 <- meta[-c(19:48),]
sapply(meta2, class)
##             Site              ARC           Season     Sample_dates 
##         "factor"         "factor"         "factor"         "factor" 
##      Pea_variety             Plot    season_precip        irrgation 
##         "factor"         "factor"        "numeric"        "numeric" 
## total_precip_irr     sample_depth             Date          Tillage 
##        "numeric"         "factor"         "factor"         "factor" 
##        prev_crop      grain_yield        elevation              lat 
##         "factor"        "numeric"         "factor"         "factor" 
##              lon   Organic_Matter Moisture_Content  Nitrate_Nitrite 
##         "factor"        "numeric"        "numeric"        "numeric" 
##          Ammonia    Av_Phosphorus     Av_Potassium   Sulfate_Sulfur 
##        "numeric"        "numeric"        "numeric"        "numeric" 
##               pH            Boron          Arsenic           Barium 
##        "numeric"        "numeric"        "numeric"        "numeric" 
##          Cadmium          Calcium         Chromium           Cobalt 
##        "numeric"        "numeric"        "numeric"        "numeric" 
##           Copper             Iron             Lead        Magnesium 
##        "numeric"        "numeric"        "numeric"        "numeric" 
##        Manganese       Molybdenum           Nickel       Phosphorus 
##        "numeric"        "numeric"        "numeric"        "numeric" 
##        Potassium           Sodium           Sulfur         Vanadium 
##        "numeric"        "numeric"        "numeric"        "numeric" 
##             Zinc 
##        "numeric"

Convert to matrix

OTU_nifH_m<-as.matrix(OTU_nifH)
tax_nifH_m<-as.matrix(tax_nifH)
meta_m<-as.matrix(meta2)

class(OTU_nifH_m)
## [1] "matrix"
class(tax_nifH_m)
## [1] "matrix"
class(meta_m)
## [1] "matrix"

Make phyloseq object

OTUnifH = otu_table(OTU_nifH_m, taxa_are_rows = TRUE)
TAXnifH = tax_table(tax_nifH_m)

physeq_nifH = phyloseq(OTUnifH, TAXnifH)

Get physeq info

physeq_nifH
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 8821 taxa and 101 samples ]
## tax_table()   Taxonomy Table:    [ 8821 taxa by 8 taxonomic ranks ]

Add meta data to phyloseq object

meta_phy <- sample_data(meta2)
sample_names(meta_phy)
##  [1] "JZ032" "JZ031" "JZ030" "JZ034" "JZ033" "JZ035" "JZ040" "JZ046"
##  [9] "JZ042" "JZ044" "JZ038" "JZ036" "JZ041" "JZ047" "JZ043" "JZ045"
## [17] "JZ039" "JZ037" "JZ081" "JZ078" "JZ082" "JZ079" "JZ080" "JZ083"
## [25] "JZ105" "JZ107" "JZ103" "JZ102" "JZ106" "JZ104" "JZ084" "JZ085"
## [33] "JZ086" "JZ087" "JZ088" "JZ089" "JZ090" "JZ091" "JZ092" "JZ093"
## [41] "JZ094" "JZ095" "JZ096" "JZ097" "JZ098" "JZ099" "JZ100" "JZ101"
## [49] "JZ112" "JZ109" "JZ110" "JZ111" "JZ108" "JZ113"
physeq_nifH<-merge_phyloseq(physeq_nifH, meta_phy)
physeq_nifH
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 8821 taxa and 54 samples ]
## sample_data() Sample Data:       [ 54 samples by 45 sample variables ]
## tax_table()   Taxonomy Table:    [ 8821 taxa by 8 taxonomic ranks ]

Trim data

Trim data to exclude OTUs that are not in any samples

Source of trim protocol http://evomics.org/wp-content/uploads/2016/01/phyloseq-Lab-01-Answers.html#taxa-total-counts-histogram

tdt_nifH = data.table(tax_table(physeq_nifH),
                 TotalCounts = taxa_sums(physeq_nifH),
                 OTU = taxa_names(physeq_nifH))
ggplot(tdt_nifH, aes(TotalCounts)) + 
  geom_histogram() + 
  ggtitle("Histogram of Total Counts nifH")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# How many OTUS have low count (Rare)?
tdt_nifH[(TotalCounts <= 0), .N]#zero count
## [1] 1691
tdt_nifH[(TotalCounts <= 1), .N]#single count
## [1] 2337
tdt_nifH[(TotalCounts <= 2), .N]#double count
## [1] 2778
# taxa cumulative sum
taxcumsum_nifH = tdt_nifH[, .N, by = TotalCounts]
setkey(taxcumsum_nifH, TotalCounts)
taxcumsum_nifH[, CumSum := cumsum(N)]
# Define the plot
pCumSum_nifH = ggplot(taxcumsum_nifH, aes(TotalCounts, CumSum)) + 
  geom_point() +
  xlab("Filtering Threshold, Minimum Total Counts") +
  ylab("OTUs Filtered") +
  ggtitle("OTUs that would be filtered vs. the minimum count threshold (nifH)")
pCumSum_nifH

Zoom in

pCumSum_nifH + xlim(0, 100)
## Warning: Removed 856 rows containing missing values (geom_point).

mdt_nifH = fast_melt(physeq_nifH)
prevdtnifH = mdt_nifH[, list(Prevalence = sum(count > 0), 
                    TotalCounts = sum(count)),
             by = TaxaID]
ggplot(prevdtnifH, aes(Prevalence)) + 
  geom_histogram() + 
  ggtitle("Histogram of Taxa Prevalence nifH")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# How many OTUS have low prevelance (Rare)?
prevdtnifH[(Prevalence <= 0), .N]#zero 
## [1] 1691
prevdtnifH[(Prevalence <= 1), .N]#single 
## [1] 3277
prevdtnifH[(Prevalence <= 2), .N]#double 
## [1] 4476
prevcumsumnifH = prevdtnifH[, .N, by = Prevalence]
setkey(prevcumsumnifH, Prevalence)
prevcumsumnifH[, CumSum := cumsum(N)]
pPrevCumSumnifH = ggplot(prevcumsumnifH, aes(Prevalence, CumSum)) + 
  geom_point() +
  xlab("Filtering Threshold, Prevalence") +
  ylab("OTUs Filtered") +
  ggtitle("OTUs that would be filtered vs. the minimum count threshold")
pPrevCumSumnifH

ggplot(prevdtnifH, aes(Prevalence, TotalCounts)) + 
  geom_point(size = 2, alpha = 0.5) + 
  scale_y_log10()

*** ###Trimming ####Remove less than doublets in data and prevlant in 5% of the sample ####First transform to realtive abundance

physeq_nifH_trim = filter_taxa(physeq_nifH, function(x) sum(x > 2) > (0.05*length(x)),
                               TRUE)
physeq_nifH_trim
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2107 taxa and 54 samples ]
## sample_data() Sample Data:       [ 54 samples by 45 sample variables ]
## tax_table()   Taxonomy Table:    [ 2107 taxa by 8 taxonomic ranks ]

To simplfy ordination and save time we will trim the OTUs more

Remove OTUs that do not show appear more than 5 times in more than 10th of the samples

wh0_n = genefilter_sample(physeq_nifH_trim, filterfun_sample(function(x) x > 5), 
                          A=0.1*nsamples(physeq_nifH_trim))
physeq_nifH_ord = prune_taxa(wh0_n, physeq_nifH_trim)
physeq_nifH_ord
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 804 taxa and 54 samples ]
## sample_data() Sample Data:       [ 54 samples by 45 sample variables ]
## tax_table()   Taxonomy Table:    [ 804 taxa by 8 taxonomic ranks ]

Transform to even sampling depth

physeq_nifH_ord = transform_sample_counts(physeq_nifH_ord, function(x) 1E6 * x/sum(x))
physeq_nifH_ord
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 804 taxa and 54 samples ]
## sample_data() Sample Data:       [ 54 samples by 45 sample variables ]
## tax_table()   Taxonomy Table:    [ 804 taxa by 8 taxonomic ranks ]

We have removed the majorty of the low abundance data with a remaining 804 taxa which make the data analysis much more managable.



Alpha Diveristy Analysis

Bar plots

physeq_nifH_ord_2_phylum <- tax_glom(physeq_nifH_ord, "phylum")
plot_bar(physeq_nifH_ord_2_phylum, x = "Site", fill = "phylum")+
geom_bar(aes(fill = phylum), color = "black", stat = "identity", position = "stack", 
         show.legend = TRUE)+
scale_fill_manual(name = "Phylum", 
                   values=c("#7fc97f", "#beaed4", "#fdc086", "#ffff99", "#386cb0", 
                            "#f0027f",'#e31a1c'),
                   #breaks=c("p_actinobacteria", "p_bacteroidetes", "p_cyanobacteria",
                         #"p_firmicutes", "p_proteobacteria", "p_verrucomicrobia"),
                   labels = c("Actinobacteria", "Bacteroidetes", "Cyanobacteria", 
                              "Firmicutes", "Proteobacteria", "Verrucomicrobia"),
                  guide = guide_legend(reverse = TRUE)
                    )+
ggtitle("nifH Phylum Relative Abundance by Site")+
ylab("Relative Abundance (%)")+
theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_blank(),
      panel.background = element_blank())

physeq_nifH_ord_2_phylum <- tax_glom(physeq_nifH_ord, "phylum")
plot_bar(physeq_nifH_ord_2_phylum, x = "Site", fill = "phylum")+
geom_bar(aes(fill = phylum), stat = "identity", position = "stack", show.legend = TRUE)+
scale_fill_manual(name = "Phylum", 
                   values=c("#7fc97f", "#beaed4", "#fdc086", "#ffff99", 
                            "#386cb0", "#f0027f"),
                   #breaks=c("p_actinobacteria", "p_bacteroidetes", "p_cyanobacteria",
                         #"p_firmicutes", "p_proteobacteria", "p_verrucomicrobia"),
                   labels = c("Actinobacteria", "Bacteroidetes", "Cyanobacteria",
                              "Firmicutes", "Proteobacteria", "Verrucomicrobia"),
                  guide = guide_legend(reverse = TRUE)
                    )+
ggtitle("nifH Phylum Relative Abundance by Site")+
ylab("Relative Abundance (%)")+
theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_blank(), 
      panel.background = element_blank())

Reduce phylums that contribute the less than 5% of total abundance into one group

physeq_nifH_ord_1 = transform_sample_counts(physeq_nifH_ord, function(x) x / sum(x) )
physeq_nifH_ord_phylum <- tax_glom(physeq_nifH_ord_1, "phylum")
data_nifH_phylum <- psmelt(physeq_nifH_ord_phylum)
data_nifH_phylum$phylum<-as.character(data_nifH_phylum$phylum)
data_nifH_phylum$phylum[data_nifH_phylum$Abundance<0.05]<-"<5% abdund"

#list new phylums 
unique(data_nifH_phylum$phylum)
## [1] "p__proteobacteria" "p__firmicutes"     "p__actinobacteria"
## [4] "<5% abdund"
ggplot(data = data_nifH_phylum, aes(x = Site, y = Abundance, fill = phylum))+
geom_bar(aes(fill = phylum), stat = "identity", position = "stack", show.legend = TRUE)+
scale_fill_manual(name = "Phylum", 
                   values=c( "#33a02c",  "#b2df8a","#a6cee3", "#1f78b4"),
                   breaks=c("p__proteobacteria", "p__firmicutes", 
                            "p__actinobacteria","<5% abdund" ),
                   labels = c( "Proteobacteria", "Firmicutes", 
                               "Actinobacteria", "<5% Abundance"),
                  guide = guide_legend(reverse = TRUE)
                    )+
ggtitle("nifH Phylum Relative Abundance by Site")+
ylab("Relative Abundance (%)")+
scale_x_discrete(labels = c("Conrad", "Corvallis", "Huntley Dryland", 
                            "Huntley Irrigated", "Kalispell", "Mocassin",
                            "Sidney Dryland", "Sidney Irrigated", "Richland"))+
theme(axis.text.x = element_text(angle = 90, hjust = 1), 
      axis.text.y = element_blank(), panel.background = element_blank())

Publish plot to tiff

## png 
##   2

Different Colors

ggplot(data = data_nifH_phylum, aes(x = Site, y = Abundance, fill = phylum))+
geom_bar(aes(fill = phylum), stat = "identity", position = "stack", show.legend = TRUE)+
scale_fill_manual(name = "Phylum", 
                   values=c( "#33a02c",  "#b2df8a","#a6cee3", "#ff7f00"),
                   breaks=c("p__proteobacteria", "p__firmicutes", 
                            "p__actinobacteria","<5% abdund" ),
                   labels = c( "Proteobacteria", "Firmicutes", 
                               "Actinobacteria", "<5% Abundance"),
                  guide = guide_legend(reverse = TRUE)
                    )+
ggtitle("nifH Phylum Relative Abundance by Site")+
ylab("Relative Abundance")+
scale_x_discrete(labels = c("Conrad", "Corvallis", 
                            "Huntley Dryland", "Huntley Irrigated", "Kalispell", 
                            "Mocassin", "Sidney Dryland", "Sidney Irrigated", "Richland"))+
theme(axis.text.x = element_text(angle = 90, hjust = 1), 
      axis.text.y = element_blank(), panel.background = element_blank())

nifH genus barplot

physeq_nifH_ord_1 = transform_sample_counts(physeq_nifH_ord, function(x) x / sum(x) )
physeq_nifH_ord_genus <- tax_glom(physeq_nifH_ord_1, "genus")
data_nifH_genus <- psmelt(physeq_nifH_ord_genus)
data_nifH_genus$genus<-as.character(data_nifH_genus$genus)
data_nifH_genus$genus[data_nifH_genus$Abundance<0.1]<-"<10% abdund"
unique(data_nifH_genus$genus)
##  [1] "g__paenibacillus"     "g__rhizobium"         "g__geobacter"        
##  [4] "g__rhodopseudomonas"  "g__bradyrhizobium"    "g__gluconacetobacter"
##  [7] "g__dechloromonas"     "g__agrobacterium"     "g__ideonella"        
## [10] "g__azospirillum"      "g__arthrobacter"      "<10% abdund"
ggplot(data = data_nifH_genus, aes(x = Site, y = Abundance, fill = genus))+
geom_bar(aes(fill = genus), stat = "identity", position = "stack", show.legend = TRUE)+
scale_fill_manual(name = "Genus", 
                  values=c('#fdbf6f','#a6cee3','#b2df8a','#fb9a99','#e31a1c','#ff7f00',                                               '#cab2d6','#6a3d9a','#ffff99','#80b1d3','#1f78b4', '#33a02c'),
                  breaks=c("g__paenibacillus", "g__rhizobium", "g__geobacter", 
                           "g__rhodopseudomonas","g__bradyrhizobium",
                           "g__gluconacetobacter","g__dechloromonas", 
                           "g__agrobacterium", "g__ideonella",  
                           "g__azospirillum","g__arthrobacter", "<10% abdund"),
                  labels=c("Paenibacillus", "Rhizobium", "Geobacter", "Rhodopseudomonas",                                             "Bradyrhizobium", "Gluconacetobacter", "Dechloromonas",                                                    "Agrobacterium", "Ideonella", "Azospirillum","Arthrobacter", 
                            "<10% abdund"),
                  guide = guide_legend(reverse = TRUE)
                    )+
ggtitle("nifH Genus Relative Abundance by Site")+
ylab("Relative Abundance")+
scale_x_discrete(labels = c("Conrad", "Corvallis", "Huntley Dryland", 
                            "Huntley Irrigated", "Kalispell", "Mocassin", 
                            "Sidney Dryland", "Sidney Irrigated", "Richland"))+
theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_blank(),
      panel.background = element_blank())

Publish to a tiff image

## png 
##   2

Alpha diversity metrics

Use phyloseq internal packages to calculate the alpha diverisity

plot_richness(physeq_nifH_trim)

Simplify to just Observed and Chao1 and Shannon

plot_richness(physeq_nifH_trim, measures = c("Observed","Chao1", "Shannon"), color = "Site")

rich_nifH<-estimate_richness(physeq_nifH_trim, split = TRUE)


#First merge data sets with meta2 
meta2$sample_names<-rownames(meta2)
rich_nifH$sample_names<-rownames(rich_nifH)
meta_nifH<-merge(meta2, rich_nifH, by = "sample_names")
rownames(meta_nifH)<-meta_nifH$sample_names
meta_nifH<-meta_nifH[,-1]
head(meta_nifH)[,49:54]
##             ACE   se.ACE  Shannon   Simpson InvSimpson   Fisher
## JZ030 1061.8638 16.02588 3.840892 0.9370300  15.880569 130.0829
## JZ031 1015.7345 15.87722 3.786340 0.9379526  16.116700 120.5734
## JZ032  985.3458 15.23833 3.241100 0.8739992   7.936459 116.9659
## JZ033 1077.3023 16.24837 3.968293 0.9463960  18.655339 133.1658
## JZ034 1044.4015 15.97121 3.962209 0.9532236  21.378282 129.1375
## JZ035 1123.6740 16.58229 3.653639 0.9231705  13.015828 133.4697

Make boxplots with Observed Shannon chao1 and inverse simpson

#use ggpubr for plot 
nifH_Observ<-ggboxplot(meta_nifH, x = "Site", y = "Observed",
   rug = TRUE,
   fill = "Site", xlab = " ",
   palette = farm_col_paired)+
  rremove("x.text")
  
nifH_Shannon<-ggboxplot(meta_nifH, x = "Site", y = "Shannon",
    rug = TRUE,
   fill = "Site", xlab = " ",
   palette = farm_col_paired)+
  rremove("x.text")
  
nifH_Chao<- ggboxplot(meta_nifH, x = "Site", y = "Chao1",
    rug = TRUE,
   fill = "Site", xlab = " ",
   palette = farm_col_paired)+
  rremove("x.text")

nifH_InvSim<- ggboxplot(meta_nifH, x = "Site", y = "InvSimpson",
   rug = TRUE,
   fill = "Site", xlab = " ",
   palette = farm_col_paired)+
  rremove("x.text")

alpha_nifH_fig<-ggarrange(nifH_Observ, nifH_Shannon, nifH_Chao, nifH_InvSim, ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom")

annotate_figure(alpha_nifH_fig, top = text_grob("Alpha Diversity of nifH by Site", size = 20))   

Publish to .tiff

## png 
##   2

Observed- total observed OTUs

Chao1- estimate diversity and assumes that the number of observations

for a taxa has a Poisson distribution and corrects for variance

Shannon- # of OTUs (richness) scaled to the evenness

Simpson- scale of dominace probabilty of any two indviduals drawn at random beloging to the same species

Source: http://biology.kenyon.edu/courses/biol229/diversity.pdf

Just shannon plot by site

#use ggpubr for plot 
ggboxplot(meta_nifH, x = "Site", y = "Shannon",
   add = "mean", rug = TRUE,
   fill = "Site",
   title = "nifH Shannon Diversity by Site", palette = farm_col_paired, legend = "right")+
  stat_compare_means(label.y = 4.3,  p.adjust.method = "bonferroni")+
  rotate_x_text()
## Warning: Ignoring unknown parameters: p.adjust.method

Determine if distrubution is normal for each diversity test

Used the following protocol :

https://rpubs.com/dillmcfarlan/R_microbiotaSOP

#Create 2x2 plot environment so that we can see all 4 metrics at once. 
par(mfrow = c(2, 2))

#Then plot each metric.
hist(rich_nifH$Shannon, main="Shannon diversity", xlab="", breaks=10)
hist(rich_nifH$InvSimpson, main="Inverse Simpson diversity", xlab="", breaks=15)
hist(rich_nifH$Chao1, main="Chao richness", xlab="", breaks=15)
hist(rich_nifH$ACE, main="ACE richness", xlab="", breaks=15)

Test for normalcy using the shapiro test. The null hypothesis for this test is that the data are normally distributed, if the p-value is greater than 0.05, then the null hypothesis is not rejected.

shapiro.test(rich_nifH$Shannon)
## 
##  Shapiro-Wilk normality test
## 
## data:  rich_nifH$Shannon
## W = 0.97136, p-value = 0.2211
shapiro.test(rich_nifH$InvSimpson)
## 
##  Shapiro-Wilk normality test
## 
## data:  rich_nifH$InvSimpson
## W = 0.8798, p-value = 6.121e-05
shapiro.test(rich_nifH$Chao1)
## 
##  Shapiro-Wilk normality test
## 
## data:  rich_nifH$Chao1
## W = 0.94469, p-value = 0.01469
shapiro.test(rich_nifH$ACE)
## 
##  Shapiro-Wilk normality test
## 
## data:  rich_nifH$ACE
## W = 0.93376, p-value = 0.005166

Use ANOVA on alpha diveristy metrics for Shannon variables becuase the other are not normal

Shannon first

aov_shannon_site_nifH <- aov(Shannon ~ Site, meta_nifH)
summary(aov_shannon_site_nifH)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Site         8 15.624  1.9531   15.43 1.28e-10 ***
## Residuals   45  5.696  0.1266                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Site location has a signifcant effect on shannon diversity

Correct for multiple comparisons

TukeyHSD(aov_shannon_site_nifH)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Shannon ~ Site, data = meta_nifH)
## 
## $Site
##                                           diff          lwr          upr
## Corvallis-Conrad                    0.29300868 -0.376019931  0.962037285
## Huntley Dryland -Conrad            -0.50522512 -1.174253731  0.163803484
## Huntley Irrigated-Conrad            0.86418474  0.195156134  1.533213349
## Kalispell-Conrad                    0.66189448 -0.007134125  1.330923091
## Mocassin-Conrad                    -0.44378350 -1.112812113  0.225245103
## Richland-Conrad                    -0.68161691 -1.350645518 -0.012588302
## Sidney Dryland-Conrad              -0.66014194 -1.329170551  0.008886665
## Sidney Irrigated-Conrad             0.03893720 -0.630091409  0.707965807
## Huntley Dryland -Corvallis         -0.79823380 -1.467262408 -0.129205192
## Huntley Irrigated-Corvallis         0.57117606 -0.097852543  1.240204673
## Kalispell-Corvallis                 0.36888581 -0.300142802  1.037914414
## Mocassin-Corvallis                 -0.73679218 -1.405820790 -0.067763574
## Richland-Corvallis                 -0.97462559 -1.643654195 -0.305596979
## Sidney Dryland-Corvallis           -0.95315062 -1.622179227 -0.284122012
## Sidney Irrigated-Corvallis         -0.25407148 -0.923100086  0.414957130
## Huntley Irrigated-Huntley Dryland   1.36940986  0.700381257  2.038438473
## Kalispell-Huntley Dryland           1.16711961  0.498090998  1.836148214
## Mocassin-Huntley Dryland            0.06144162 -0.607586989  0.730470227
## Richland-Huntley Dryland           -0.17639179 -0.845420395  0.492636821
## Sidney Dryland-Huntley Dryland     -0.15491682 -0.823945427  0.514111789
## Sidney Irrigated-Huntley Dryland    0.54416232 -0.124866286  1.213190930
## Kalispell-Huntley Irrigated        -0.20229026 -0.871318867  0.466738349
## Mocassin-Huntley Irrigated         -1.30796825 -1.976996854 -0.638939638
## Richland-Huntley Irrigated         -1.54580165 -2.214830260 -0.876773044
## Sidney Dryland-Huntley Irrigated   -1.52432668 -2.193355292 -0.855298076
## Sidney Irrigated-Huntley Irrigated -0.82524754 -1.494276151 -0.156218935
## Mocassin-Kalispell                 -1.10567799 -1.774706596 -0.436649380
## Richland-Kalispell                 -1.34351139 -2.012540001 -0.674482785
## Sidney Dryland-Kalispell           -1.32203643 -1.991065033 -0.653007818
## Sidney Irrigated-Kalispell         -0.62295728 -1.291985892  0.046071324
## Richland-Mocassin                  -0.23783341 -0.906862013  0.431195203
## Sidney Dryland-Mocassin            -0.21635844 -0.885387046  0.452670170
## Sidney Irrigated-Mocassin           0.48272070 -0.186307904  1.151749312
## Sidney Dryland-Richland             0.02147497 -0.647553640  0.690503576
## Sidney Irrigated-Richland           0.72055411  0.051525501  1.389582717
## Sidney Irrigated-Sidney Dryland     0.69907914  0.030050533  1.368107749
##                                        p adj
## Corvallis-Conrad                   0.8815206
## Huntley Dryland -Conrad            0.2786584
## Huntley Irrigated-Conrad           0.0035559
## Kalispell-Conrad                   0.0545026
## Mocassin-Conrad                    0.4476016
## Richland-Conrad                    0.0428516
## Sidney Dryland-Conrad              0.0556618
## Sidney Irrigated-Conrad            0.9999999
## Huntley Dryland -Corvallis         0.0091627
## Huntley Irrigated-Corvallis        0.1498601
## Kalispell-Corvallis                0.6844416
## Mocassin-Corvallis                 0.0211504
## Richland-Corvallis                 0.0006683
## Sidney Dryland-Corvallis           0.0009315
## Sidney Irrigated-Corvallis         0.9435800
## Huntley Irrigated-Huntley Dryland  0.0000011
## Kalispell-Huntley Dryland          0.0000308
## Mocassin-Huntley Dryland           0.9999977
## Richland-Huntley Dryland           0.9940451
## Sidney Dryland-Huntley Dryland     0.9975503
## Sidney Irrigated-Huntley Dryland   0.1957237
## Kalispell-Huntley Irrigated        0.9854512
## Mocassin-Huntley Irrigated         0.0000030
## Richland-Huntley Irrigated         0.0000001
## Sidney Dryland-Huntley Irrigated   0.0000001
## Sidney Irrigated-Huntley Irrigated 0.0062523
## Mocassin-Kalispell                 0.0000835
## Richland-Kalispell                 0.0000017
## Sidney Dryland-Kalispell           0.0000024
## Sidney Irrigated-Kalispell         0.0858435
## Richland-Mocassin                  0.9611572
## Sidney Dryland-Mocassin            0.9778887
## Sidney Irrigated-Mocassin          0.3356109
## Sidney Dryland-Richland            1.0000000
## Sidney Irrigated-Richland          0.0261606
## Sidney Irrigated-Sidney Dryland    0.0344455

Plot Irrigation managment against diffrent diveristies

nifH_irr_Observ<-ggboxplot(meta_nifH, x = "Plot", y = "Observed",
   rug = TRUE,
   fill = "Plot", xlab = " ",
   palette = farm_col_paired)+
  rremove("x.text")
  
nifH_irr_Shannon<-ggboxplot(meta_nifH, x = "Plot", y = "Shannon",
   rug = TRUE,
   fill = "Plot", xlab = " ",
   palette = farm_col_paired)+
  rremove("x.text")
  
nifH_irr_Chao<- ggboxplot(meta_nifH, x = "Plot", y = "Chao1",
    rug = TRUE,
   fill = "Plot", xlab = " ",
   palette = farm_col_paired)+
  rremove("x.text")

nifH_irr_InvSim<- ggboxplot(meta_nifH, x = "Plot", y = "Simpson",
   rug = TRUE,
   fill = "Plot", xlab = " ",
   palette = farm_col_paired)+
  rremove("x.text")

alpha_nifH_irr_fig<-ggarrange(nifH_irr_Observ, nifH_irr_Shannon, nifH_irr_Chao, nifH_irr_InvSim, ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom")

annotate_figure(alpha_nifH_irr_fig, top = text_grob("Alpha Diversity of nifH by Irrigation", size = 20))   

Looks like there is significance in the other diverisities to for irrgation (Dryland vs irrigated)

aov_observed_irr_nifH <- aov(Observed ~ Plot, meta_nifH)
summary(aov_observed_irr_nifH)
##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## Plot         1 1351360 1351360   81.29 3.28e-12 ***
## Residuals   52  864472   16624                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Signficance in observed

aov_observed_irr_nifH <- aov(Observed ~ Plot/Site, meta_nifH)
summary(aov_observed_irr_nifH)
##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## Plot         1 1351360 1351360  337.86  < 2e-16 ***
## Plot:Site    7  684482   97783   24.45 2.32e-13 ***
## Residuals   45  179990    4000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot is significant but the interaction of Plot:Site is also signficant so we cannot say there is actuall correlation

aov_Chao1_irr_nifH <- aov(Chao1 ~ Plot, meta_nifH)
summary(aov_Chao1_irr_nifH)
##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## Plot         1 1328751 1328751   62.35 1.86e-10 ***
## Residuals   52 1108257   21313                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Signifcant in chao1

aov_Chao1_irr_nifH <- aov(Chao1 ~ Plot/Site, meta_nifH)
summary(aov_Chao1_irr_nifH)
##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## Plot         1 1328751 1328751  234.21  < 2e-16 ***
## Plot:Site    7  852953  121850   21.48 2.11e-12 ***
## Residuals   45  255304    5673                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
aov_Simpson_irr_nifH <- aov(Simpson ~ Plot, meta_nifH)
summary(aov_Simpson_irr_nifH)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Plot         1 0.1043 0.10431   19.64 4.85e-05 ***
## Residuals   52 0.2762 0.00531                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
aov_Simpson_irr_nifH <- aov(Simpson ~ Plot/Site, meta_nifH)
summary(aov_Simpson_irr_nifH)
##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## Plot         1 0.10431 0.10431  23.374 1.59e-05 ***
## Plot:Site    7 0.07539 0.01077   2.413   0.0347 *  
## Residuals   45 0.20083 0.00446                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

All alpha diversities are signficant with irrgation method but they are nested in Site so we must look within site (Huntley and Sidney) later. It is hard to prove anything in the statewide survey with univariate ANOVAs

nifH_irr_Observ_scatter<-ggscatter(meta_nifH, x = "total_precip_irr", y = "Observed",
   rug = TRUE,
   add = "reg.line",
   color = "Plot", xlab = " ",
   palette = farm_col_paired)+
  rremove("x.text")
  
nifH_irr_Shannon_scatter<-ggscatter(meta_nifH, x = "total_precip_irr", y = "Shannon",
   rug = TRUE,
   add = "reg.line",
   color = "Plot", xlab = " ",
   palette = farm_col_paired)+
  rremove("x.text")
  
nifH_irr_Chao_scatter<- ggscatter(meta_nifH, x = "total_precip_irr", y = "Chao1",
   rug = TRUE,
   add = "reg.line",
   color = "Plot", xlab = " ",
   palette = farm_col_paired)+
  rremove("x.text")

nifH_irr_InvSim_scatter<- ggscatter(meta_nifH, x = "total_precip_irr", y = "Simpson",
   rug = TRUE,
   add = "reg.line",
   color = "Plot", xlab = " ",
   palette = farm_col_paired)+
  rremove("x.text")

alpha_nifH_irr_fig<-ggarrange(nifH_irr_Observ_scatter, nifH_irr_Shannon_scatter, 
                              nifH_irr_Chao_scatter, nifH_irr_InvSim_scatter,
                              ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom")

annotate_figure(alpha_nifH_irr_fig, 
                top = text_grob("Alpha Diversity of nifH by Irrigation", size = 20))   


Plot ordination

Used the following protocol from the phyloseq tutorial:

https://joey711.github.io/phyloseq/plot_ordination-examples.html

RDA Ordination

phynifH_ord_RDA <- ordinate(physeq_nifH_ord, "RDA", "bray")
plot_ordination(physeq_nifH_ord, phynifH_ord_RDA, color = "Site", shape = "Plot")+
    geom_point(size = 3)+
  scale_color_manual(values =  farm_col_paired)+
  stat_ellipse(type = "norm", linetype = 2, aes(color ="Plot"), show.legend = TRUE) +
  ggtitle("nifH RDA Ordination by Site and Irrigation")+
  theme_bw()

Irrigation is influencing the ordination of the principle components (PC1 is most likely comprsied of irrigation / other farm Managment)

This is really diffrenent from the irrigation RDA plot for 16s

Publish to the tiff

## png 
##   2

Will color with other farm managment to see if anytning is intresting.

plot_ordination(physeq_nifH_ord, phynifH_ord_RDA, color = "prev_crop", shape = "Plot")+
    geom_point(size = 3)+
  scale_color_manual(values =  farm_col_paired)+
  ggtitle("nifH RDA Ordination by Previous Crop and Irrigation")+
  theme_bw()

plot_ordination(physeq_nifH_ord, phynifH_ord_RDA, color = "Tillage", shape = "Plot")+
    geom_point(size = 3)+
  scale_color_manual(values =  farm_col_paired)+
  ggtitle("nifH RDA Ordination by Previous Crop and Irrigation")+
  theme_bw()

RDA plots are much diffrenent from the 16s versions but the two RDA plots from both 16s and nifH are still overlapping

Pea variety has no correlation or ordnaiton to bacterial community bray-curtis distance

plot_ordination(physeq_nifH_ord, phynifH_ord_RDA, color = "Pea_variety", shape = "Plot")+
    geom_point(size = 3)+
  scale_color_manual(values =  farm_col_paired)+
  ggtitle("nifH RDA Ordination by Previous Crop and Irrigation")+
  theme_bw()

####RDA has alot of overlap


NMDS ordination

phynifH_ord_NMDS <- ordinate(physeq_nifH_ord, "NMDS", "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1618779 
## Run 1 stress 0.1666296 
## Run 2 stress 0.1618753 
## ... New best solution
## ... Procrustes: rmse 0.001579567  max resid 0.008369045 
## ... Similar to previous best
## Run 3 stress 0.1887229 
## Run 4 stress 0.161878 
## ... Procrustes: rmse 0.001580814  max resid 0.008390491 
## ... Similar to previous best
## Run 5 stress 0.1618754 
## ... Procrustes: rmse 4.945684e-05  max resid 0.0003091343 
## ... Similar to previous best
## Run 6 stress 0.1618781 
## ... Procrustes: rmse 0.001584558  max resid 0.008371028 
## ... Similar to previous best
## Run 7 stress 0.1657018 
## Run 8 stress 0.1618754 
## ... Procrustes: rmse 3.746306e-05  max resid 0.0002276462 
## ... Similar to previous best
## Run 9 stress 0.1887328 
## Run 10 stress 0.1618754 
## ... Procrustes: rmse 9.792843e-05  max resid 0.0005519838 
## ... Similar to previous best
## Run 11 stress 0.1618754 
## ... Procrustes: rmse 7.118828e-05  max resid 0.0004423043 
## ... Similar to previous best
## Run 12 stress 0.1618779 
## ... Procrustes: rmse 0.001579864  max resid 0.008364055 
## ... Similar to previous best
## Run 13 stress 0.1618779 
## ... Procrustes: rmse 0.00158213  max resid 0.008386758 
## ... Similar to previous best
## Run 14 stress 0.1618781 
## ... Procrustes: rmse 0.001578848  max resid 0.008382921 
## ... Similar to previous best
## Run 15 stress 0.1618779 
## ... Procrustes: rmse 0.001582478  max resid 0.008396537 
## ... Similar to previous best
## Run 16 stress 0.1618754 
## ... Procrustes: rmse 4.310258e-05  max resid 0.0002645992 
## ... Similar to previous best
## Run 17 stress 0.1618754 
## ... Procrustes: rmse 9.174953e-05  max resid 0.0004160894 
## ... Similar to previous best
## Run 18 stress 0.1618754 
## ... Procrustes: rmse 4.384431e-05  max resid 0.0002702027 
## ... Similar to previous best
## Run 19 stress 0.1618753 
## ... Procrustes: rmse 2.118693e-05  max resid 0.00010052 
## ... Similar to previous best
## Run 20 stress 0.1618754 
## ... Procrustes: rmse 2.938269e-05  max resid 0.0001766266 
## ... Similar to previous best
## *** Solution reached
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS)

phynifH_ord_NMDS
## 
## Call:
## metaMDS(comm = veganifyOTU(physeq), distance = distance) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(veganifyOTU(physeq))) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.1618753 
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'

After stress test run, we get a value of 0.07 which is considered good, anything below 0.2 is acceptable.


####Plot NMDS with Site and Irrigation

plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "Site")+
  geom_point(size = 3)+
  scale_color_manual(values =  farm_col_paired)+
  ggtitle("nifH NMDS Ordination by Site")+
  theme_bw()

#export physeq OTU for network analysis and saving
nifH_trim_OTU <- as(otu_table(physeq_nifH_ord), "matrix")
ID <- rownames(nifH_trim_OTU)
nifH_trim_OTU_1 <- cbind(ID, nifH_trim_OTU)
rownames(nifH_trim_OTU_1)<-c()

write.table(nifH_trim_OTU_1, file = "nifH_trim_OTU.txt", sep = "\t", row.names = FALSE)
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "Pea_variety")+
  geom_point(size = 3)+
  scale_color_manual(values =  farm_col_paired)+
  ggtitle("nifH NMDS Ordination by Pea Variety")+
  theme_bw()

Pea variety does not influence the micriobiome compostion in nifH as well

plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, shape = "Plot", color = "Site")+
  geom_point(size = 3)+
  scale_color_manual(values =  farm_col_paired)+
  ggtitle("nifH NMDS Ordination by Irrigation method")+
  stat_ellipse(type = "norm", linetype = 2, aes(color = "Plot")) +
  theme_bw()

Publish to tiff

## png 
##   2

Total water added

plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "total_precip_irr")+
  geom_point(size = 3)+
  scale_color_gradient(low='#05D9F6', high='#5011D1')+
  ggtitle("nifH NMDS Ordination by Total Precipitation")+
  theme_bw()

plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "Moisture_Content")+
  geom_point(size = 3)+
  scale_color_continuous(low='#05D9F6', high='#5011D1')+
  ggtitle("nifH NMDS Ordination by Mositure Content")+
  theme_bw()

Again total water does not ordinate well, how can we explain a diffrence between irrgation and rainfall?

Lets try some of the other variables like till, fallow, and previous crop.

plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS,color = "Tillage")+
  geom_point(size = 3)+
  scale_color_manual(values =  farm_col_dark, breaks=c("Conventional", "Culti-roller", "No_till"),
                    labels=c("Conventional", "Culti Roller", "No Till"))+
  ggtitle("nifH NMDS Ordination by Tillage")+
  theme_bw()

Publish to .tiff

## png 
##   2
plot_ordination(physeq_nifH_ord, phynifH_ord_NMDS, color = "prev_crop")+
  geom_point(size = 3)+
  scale_color_manual(values =  c("#1B9E77", "#D95F02", "#7570B3", "#66a61E"),
                    name = "Previous Crop",
                    breaks=c("barley",  "Chem_fallow", "Spring_wheat", "winter_wheat"),
                    labels=c("Barley", "Chemical Fallow", "Spring Wheat", "Winter Wheat"))+
  ggtitle("nifH NMDS Ordination by Previous Crop")+
  theme_bw()

Publish to .tiff

## png 
##   2

Beta dispersions

Test the diffrences in group homogeneities. Do our farm managment factors effect the homogeneitiey of the bray curtis distance?


####If a group (Site) in the MDS space are close but have diffrenent dispersion you could have a signifcant results when it is only a diffrence in dispersion.

Anderson (2006)-https://www.ncbi.nlm.nih.gov/pubmed/16542252 https://onlinelibrary.wiley.com/doi/epdf/10.1111/j.1461-0248.2006.00926.x

Irrigation beta dispersion

disp.plot <- betadisper(distance(physeq_nifH_ord, method = "bray"), meta2$Plot)
permutest(disp.plot, pairwise=TRUE, permutations=1000)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 1000
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm  Pr(>F)  
## Groups     1 0.023867 0.0238672 6.7103   1000 0.01199 *
## Residuals 52 0.184953 0.0035568                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##            dryland irrigated
## dryland                0.012
## irrigated 0.012408
boxplot(disp.plot)

plot(disp.plot, hull = FALSE, ellipse = TRUE)

disp.Till <- betadisper(distance(physeq_nifH_ord, method = "bray"), meta2$Tillage)
permutest(disp.Till, pairwise=TRUE, permutations=1000)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 1000
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm  Pr(>F)  
## Groups     2 0.039609 0.0198044 4.2834   1000 0.01898 *
## Residuals 51 0.235799 0.0046235                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##              Conventional Culti-roller No_till
## Conventional                9.9900e-04  0.0769
## Culti-roller   4.4826e-05               0.1399
## No_till        8.3245e-02   1.4099e-01
boxplot(disp.Till)

plot(disp.Till, hull = FALSE, ellipse = TRUE)

disp.prev_crop <- betadisper(distance(physeq_nifH_ord, method = "bray"), meta2$prev_crop)
permutest(disp.prev_crop, pairwise=TRUE, permutations=1000)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 1000
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq     F N.Perm   Pr(>F)    
## Groups     3 0.34683 0.115610 14.42   1000 0.000999 ***
## Residuals 50 0.40086 0.008017                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##                  barley Chem_fallow Spring_wheat winter_wheat
## barley                   8.2917e-02   9.9900e-04       0.0010
## Chem_fallow  8.1169e-02               2.9970e-03       0.0020
## Spring_wheat 3.1442e-04  6.5930e-04                    0.6963
## winter_wheat 2.7285e-05  3.6151e-05   6.7668e-01
boxplot(disp.prev_crop)

plot(disp.prev_crop, hull = FALSE, ellipse = TRUE)

TukeyHSD(disp.prev_crop)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = distances ~ group, data = df)
## 
## $group
##                                  diff        lwr         upr     p adj
## Chem_fallow-barley        -0.03352468 -0.1077209  0.04067151 0.6292204
## Spring_wheat-barley       -0.18749276 -0.2996669 -0.07531866 0.0002814
## winter_wheat-barley       -0.22719822 -0.3393723 -0.11502412 0.0000114
## Spring_wheat-Chem_fallow  -0.15396808 -0.2625802 -0.04535597 0.0023889
## winter_wheat-Chem_fallow  -0.19367354 -0.3022856 -0.08506143 0.0001047
## winter_wheat-Spring_wheat -0.03970546 -0.1770901  0.09767919 0.8684302
disp.site <- betadisper(distance(physeq_nifH_ord, method = "bray"), meta2$Site)
permutest(disp.site, pairwise=TRUE, permutations=1000)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 1000
## 
## Response: Distances
##           Df  Sum Sq Mean Sq      F N.Perm   Pr(>F)   
## Groups     8 0.35944 0.04493 3.2652   1000 0.005994 **
## Residuals 45 0.61921 0.01376                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##                       Conrad  Corvallis Huntley Dryland  Huntley Irrigated
## Conrad                       0.03496503       0.03296703        0.06093906
## Corvallis         0.03646739                  0.68131868        0.46953047
## Huntley Dryland   0.03003786 0.67097611                         0.29470529
## Huntley Irrigated 0.06780670 0.44317758       0.29284235                  
## Kalispell         0.91870043 0.00041400       0.00027429        0.00246753
## Mocassin          0.52173883 0.07999231       0.06329790        0.16335615
## Richland          0.57388386 0.00439742       0.00283358        0.02313450
## Sidney Dryland    0.99250039 0.00143581       0.00099248        0.00635063
## Sidney Irrigated  0.32390060 0.21919454       0.17927660        0.38365390
##                    Kalispell   Mocassin   Richland Sidney Dryland
## Conrad            0.92007992 0.54045954 0.55844156     0.98601399
## Corvallis         0.00099900 0.08691309 0.00899101     0.00399600
## Huntley Dryland   0.00199800 0.07492507 0.00599401     0.00399600
## Huntley Irrigated 0.00899101 0.17082917 0.02897103     0.00699301
## Kalispell                    0.31168831 0.26973027     0.86113886
## Mocassin          0.32083698            0.81818182     0.39160839
## Richland          0.27360130 0.82416350                0.37962038
## Sidney Dryland    0.86114730 0.40276122 0.39167738               
## Sidney Irrigated  0.15341808 0.67667507 0.47279724     0.20368026
##                   Sidney Irrigated
## Conrad                      0.3167
## Corvallis                   0.2078
## Huntley Dryland             0.1738
## Huntley Irrigated           0.3956
## Kalispell                   0.1409
## Mocassin                    0.6643
## Richland                    0.4446
## Sidney Dryland              0.2138
## Sidney Irrigated
boxplot(disp.site)

plot(disp.site, hull = FALSE, ellipse = TRUE)

TukeyHSD(disp.site)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = distances ~ group, data = df)
## 
## $group
##                                             diff         lwr         upr
## Corvallis-Conrad                    0.1970866601 -0.02350369 0.417677013
## Huntley Dryland -Conrad             0.2056910662 -0.01489929 0.426281420
## Huntley Irrigated-Conrad            0.1730650455 -0.04752531 0.393655399
## Kalispell-Conrad                   -0.0092615904 -0.22985194 0.211328763
## Mocassin-Conrad                     0.0683743608 -0.15221599 0.288964714
## Richland-Conrad                     0.0514247120 -0.16916564 0.272015065
## Sidney Dryland-Conrad               0.0008764728 -0.21971388 0.221466826
## Sidney Irrigated-Conrad             0.1080798206 -0.11251053 0.328670174
## Huntley Dryland -Corvallis          0.0086044061 -0.21198595 0.229194759
## Huntley Irrigated-Corvallis        -0.0240216145 -0.24461197 0.196568739
## Kalispell-Corvallis                -0.2063482505 -0.42693860 0.014242103
## Mocassin-Corvallis                 -0.1287122993 -0.34930265 0.091878054
## Richland-Corvallis                 -0.1456619481 -0.36625230 0.074928405
## Sidney Dryland-Corvallis           -0.1962101873 -0.41680054 0.024380166
## Sidney Irrigated-Corvallis         -0.0890068395 -0.30959719 0.131583514
## Huntley Irrigated-Huntley Dryland  -0.0326260207 -0.25321637 0.187964333
## Kalispell-Huntley Dryland          -0.2149526566 -0.43554301 0.005637697
## Mocassin-Huntley Dryland           -0.1373167054 -0.35790706 0.083273648
## Richland-Huntley Dryland           -0.1542663542 -0.37485671 0.066323999
## Sidney Dryland-Huntley Dryland     -0.2048145934 -0.42540495 0.015775760
## Sidney Irrigated-Huntley Dryland   -0.0976112456 -0.31820160 0.122979108
## Kalispell-Huntley Irrigated        -0.1823266359 -0.40291699 0.038263717
## Mocassin-Huntley Irrigated         -0.1046906847 -0.32528104 0.115899669
## Richland-Huntley Irrigated         -0.1216403336 -0.34223069 0.098950020
## Sidney Dryland-Huntley Irrigated   -0.1721885728 -0.39277893 0.048401781
## Sidney Irrigated-Huntley Irrigated -0.0649852250 -0.28557558 0.155605128
## Mocassin-Kalispell                  0.0776359512 -0.14295440 0.298226305
## Richland-Kalispell                  0.0606863024 -0.15990405 0.281276656
## Sidney Dryland-Kalispell            0.0101380632 -0.21045229 0.230728417
## Sidney Irrigated-Kalispell          0.1173414110 -0.10324894 0.337931764
## Richland-Mocassin                  -0.0169496488 -0.23754000 0.203640704
## Sidney Dryland-Mocassin            -0.0674978880 -0.28808824 0.153092465
## Sidney Irrigated-Mocassin           0.0397054598 -0.18088489 0.260295813
## Sidney Dryland-Richland            -0.0505482392 -0.27113859 0.170042114
## Sidney Irrigated-Richland           0.0566551086 -0.16393524 0.277245462
## Sidney Irrigated-Sidney Dryland     0.1072033478 -0.11338701 0.327793701
##                                        p adj
## Corvallis-Conrad                   0.1134000
## Huntley Dryland -Conrad            0.0849905
## Huntley Irrigated-Conrad           0.2342545
## Kalispell-Conrad                   1.0000000
## Mocassin-Conrad                    0.9829927
## Richland-Conrad                    0.9974322
## Sidney Dryland-Conrad              1.0000000
## Sidney Irrigated-Conrad            0.8019371
## Huntley Dryland -Corvallis         1.0000000
## Huntley Irrigated-Corvallis        0.9999912
## Kalispell-Corvallis                0.0830921
## Mocassin-Corvallis                 0.6169037
## Richland-Corvallis                 0.4537416
## Sidney Dryland-Corvallis           0.1166882
## Sidney Irrigated-Corvallis         0.9217819
## Huntley Irrigated-Huntley Dryland  0.9999076
## Kalispell-Huntley Dryland          0.0613855
## Mocassin-Huntley Dryland           0.5332525
## Richland-Huntley Dryland           0.3765001
## Sidney Dryland-Huntley Dryland     0.0875794
## Sidney Irrigated-Huntley Dryland   0.8754541
## Kalispell-Huntley Irrigated        0.1797308
## Mocassin-Huntley Irrigated         0.8276676
## Richland-Huntley Irrigated         0.6843270
## Sidney Dryland-Huntley Irrigated   0.2399571
## Sidney Irrigated-Huntley Irrigated 0.9876732
## Mocassin-Kalispell                 0.9633505
## Richland-Kalispell                 0.9921001
## Sidney Dryland-Kalispell           1.0000000
## Sidney Irrigated-Kalispell         0.7237181
## Richland-Mocassin                  0.9999994
## Sidney Dryland-Mocassin            0.9843152
## Sidney Irrigated-Mocassin          0.9995991
## Sidney Dryland-Richland            0.9977226
## Sidney Irrigated-Richland          0.9950077
## Sidney Irrigated-Sidney Dryland    0.8087566


####The results from the betadispersion show that we have a signficant diffrence in the heterogeneity of our sites due to each of the farm amnagament factors. We see that there are sites that do not have signifcant diffrences and others that do. Beta dispersion cannot except models so we cannot detect if the diffrences are due to nestedness (prev_crop/Site).
####Do to the design of the experiment it will be hard to determine if the farm managment practices are responsible for the variation in bacterial population.


PERMANOVA (adoins)


####We will test the significance of the farm managment using the permutated-mulitvariate-ANOVA function in vegan called adonis. Adonis can test models and nestedness though the above beta dispersion test showed that most of the signficance is due to the diseprsion we cannot say much. But we can show we have the location effect.

adonis(distance(physeq_nifH_ord, method = "bray") ~Plot*prev_crop*Tillage*Pea_variety, data = meta2, permutations = 1000)
## 
## Call:
## adonis(formula = distance(physeq_nifH_ord, method = "bray") ~      Plot * prev_crop * Tillage * Pea_variety, data = meta2, permutations = 1000) 
## 
## Permutation: free
## Number of permutations: 1000
## 
## Terms added sequentially (first to last)
## 
##                          Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)    
## Plot                      1    1.4006 1.40059  3.8638 0.07104 0.000999 ***
## prev_crop                 2    2.5341 1.26705  3.4954 0.12853 0.000999 ***
## Tillage                   2    3.0720 1.53601  4.2374 0.15581 0.000999 ***
## Pea_variety               6    1.3480 0.22466  0.6198 0.06837 1.000000    
## Plot:Tillage              1    1.5247 1.52471  4.2062 0.07733 0.000999 ***
## Plot:Pea_variety          6    0.9169 0.15282  0.4216 0.04651 1.000000    
## prev_crop:Pea_variety     9    1.7119 0.19021  0.5247 0.08683 1.000000    
## Tillage:Pea_variety      10    2.1740 0.21740  0.5997 0.11026 1.000000    
## Plot:Tillage:Pea_variety  4    0.6841 0.17103  0.4718 0.03470 0.999001    
## Residuals                12    4.3499 0.36249         0.22063             
## Total                    53   19.7163                 1.00000             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(distance(physeq_nifH_ord, method = "bray") ~(Plot*prev_crop*Tillage*Pea_variety)/Site, data = meta2, permutations = 1000)
## 
## Call:
## adonis(formula = distance(physeq_nifH_ord, method = "bray") ~      (Plot * prev_crop * Tillage * Pea_variety)/Site, data = meta2,      permutations = 1000) 
## 
## Permutation: free
## Number of permutations: 1000
## 
## Terms added sequentially (first to last)
## 
##                                         Df SumsOfSqs MeanSqs F.Model
## Plot                                     1    1.4006       1       0
## prev_crop                                2    2.5341       1       0
## Tillage                                  2    3.0720       2       0
## Pea_variety                              6    1.3480       0       0
## Plot:Tillage                             1    1.5247       2       0
## Plot:Pea_variety                         6    0.9169       0       0
## prev_crop:Pea_variety                    9    1.7119       0       0
## Tillage:Pea_variety                     10    2.1740       0       0
## Plot:Tillage:Pea_variety                 4    0.6841       0       0
## Plot:prev_crop:Tillage:Pea_variety:Site 12    4.3499       0       0
## Residuals                                0    0.0000    -Inf        
## Total                                   53   19.7163                
##                                              R2 Pr(>F)
## Plot                                    0.07104      1
## prev_crop                               0.12853      1
## Tillage                                 0.15581      1
## Pea_variety                             0.06837      1
## Plot:Tillage                            0.07733      1
## Plot:Pea_variety                        0.04651      1
## prev_crop:Pea_variety                   0.08683      1
## Tillage:Pea_variety                     0.11026      1
## Plot:Tillage:Pea_variety                0.03470      1
## Plot:prev_crop:Tillage:Pea_variety:Site 0.22063      1
## Residuals                               0.00000       
## Total                                   1.00000

Data is nested within Site (Location effect) so the signficance in the bray-curtis dissimilarity with respect to plot is not significant with the data nested due to the lack of reproduction of conditions at each plot.


####Since there is issues with doing permutated anova over multivatrite data lets try to fit the chemical and farm managment data the NMDS orientation space using the envfit function in vegan


Model Selection


ENVFIT

Envfit does not like single variable values so we remove them

meta3<-meta2[,-c(3,4,10,11,27,29,35,38,42,46)]

Will remove Site catagory liek elevation, lat, long etc that do not diffrentate between site we can call these all geographical factors as they do not change between sites.

meta3<-meta3[,-c(2,11:13)]

Model fitting will be biased by chemical outliers that are in some plots the best way to avoid this is to determine the outliers (See chemical_analysis.Rmd) and remove the whole variable since functions ENVFIT and BIOENV will remove it if there are any n/a values.


####Removing Sulfate_Sulfur, Boron, Molybdenum, Potassium, Vanadium, Chromium and Sodium (From chemical_analysis.rmd)

meta3<-meta3[,-c(16,18,21,29,31)]
envfitnifH <- envfit(phynifH_ord_NMDS , meta3, na.rm = TRUE, permu= 999) 
envfitnifH
## 
## ***VECTORS
## 
##                     NMDS1    NMDS2     r2 Pr(>r)    
## season_precip     0.30096 -0.95364 0.1693  0.008 ** 
## irrgation         0.54023  0.84152 0.1767  0.007 ** 
## total_precip_irr  0.79985 -0.60021 0.1800  0.008 ** 
## grain_yield      -0.10337  0.99464 0.0052  0.870    
## Organic_Matter   -0.58271 -0.81268 0.5434  0.001 ***
## Moisture_Content  0.21634 -0.97632 0.2847  0.002 ** 
## Nitrate_Nitrite   0.01586  0.99987 0.0848  0.092 .  
## Ammonia          -0.16640 -0.98606 0.0684  0.163    
## Av_Phosphorus     0.77470  0.63232 0.2592  0.001 ***
## Av_Potassium      0.98960 -0.14383 0.3455  0.001 ***
## pH               -0.02160 -0.99977 0.0775  0.123    
## Barium           -0.09002 -0.99594 0.3593  0.001 ***
## Calcium           0.29225 -0.95634 0.2220  0.003 ** 
## Cobalt            0.00453 -0.99999 0.0831  0.115    
## Copper            0.31097 -0.95042 0.1409  0.022 *  
## Iron              0.25017 -0.96820 0.1241  0.031 *  
## Magnesium         0.22694 -0.97391 0.0394  0.342    
## Manganese         0.18760 -0.98225 0.1245  0.045 *  
## Nickel            0.63393 -0.77339 0.1124  0.033 *  
## Phosphorus        0.09644 -0.99534 0.2283  0.004 ** 
## Sulfur            0.10985 -0.99395 0.3851  0.001 ***
## Zinc              0.44396 -0.89605 0.2825  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                           NMDS1   NMDS2
## SiteConrad               0.9037  0.2518
## SiteCorvallis            0.0861  0.7313
## SiteHuntley Dryland     -0.2456  0.0338
## SiteHuntley Irrigated   -0.1426 -0.0653
## SiteKalispell           -1.6148 -0.2105
## SiteMocassin             0.0258 -0.7600
## SiteRichland             0.2990 -0.0971
## SiteSidney Dryland      -0.1793  0.4292
## SiteSidney Irrigated     0.8678 -0.3131
## Pea_varietyAC Earlystar -0.1321  0.1464
## Pea_varietyCDC Saffron  -0.0296 -0.0669
## Pea_varietyCDC Saffron   0.4876  0.1231
## Pea_varietyDelta         0.1280 -0.0638
## Pea_varietyDS Admiral    0.0304 -0.0922
## Pea_varietyMajoret      -0.0978  0.0629
## Pea_varietyNavarro      -0.0138 -0.0287
## Plotdryland              0.1607 -0.0285
## Plotirrigated           -0.2009  0.0356
## TillageConventional     -0.3088 -0.0315
## TillageCulti-roller      0.0861  0.7313
## TillageNo_till           0.1681 -0.1274
## prev_cropbarley         -0.5571  0.1518
## prev_cropChem_fallow     0.1944  0.1544
## prev_cropSpring_wheat    0.8678 -0.3131
## prev_cropwinter_wheat    0.0258 -0.7600
## 
## Goodness of fit:
##                 r2 Pr(>r)    
## Site        0.8843  0.001 ***
## Pea_variety 0.0323  0.991    
## Plot        0.0449  0.107    
## Tillage     0.1578  0.005 ** 
## prev_crop   0.4008  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

The envfit function allows us to see the correlation of our envriomental vectors to the bray-curtis species dissmilarity matrix in NMDS space. This is a loose corelation to real linear corelation but it can tell us how the NMDS orinetaion is being driven.

Try a quick plot with base r and vegan for the vectors

plot(phynifH_ord_NMDS, display = "sites")
plot(envfitnifH, p.max = 0.001 )

Organic matter is resposible for 54% of the variance

Site accounts for 88% of the variance

Season precip > Irrigation


Model Selection

Bioenv

Info:

BIOENV

“The function calculates a community dissimilarity matrix using vegdist. Then it selects all possible subsets of environmental variables, scales the variables, and calculates Euclidean distances for this subset using dist. The function finds the correlation between community dissimilarities and environmental distances, and for each size of subsets, saves the best result. There are 2^p-1 subsets of p variables, and an exhaustive search may take a very, very, very long time (parameter upto offers a partial relief).”

Mantel test

“Mantel statistic is simply a correlation between entries of two dissimilarity matrices (some use cross products, but these are linearly related). However, the significance cannot be directly assessed, because there are N(N-1)/2 entries for just N observations. Mantel developed asymptotic test, but here we use permutations of N rows and columns of dissimilarity matrix.”


OTU_nifH_trim<- as(otu_table(physeq_nifH_ord ), "matrix")
#Transpose the data to have sample names on rows
abund_tablenifH<-t(OTU_nifH_trim)

nrow(abund_tablenifH)
## [1] 54
setdiff(rownames(meta3), rownames(abund_tablenifH))
## character(0)

Our meta data and sample data match with 0 diffrence in rownames

Will use parallel processing to speed up calculations

#First detect amount of cores avalible 
detectCores()
## [1] 12
#get bray-curtis distance
abund_distnifH<-vegdist(abund_tablenifH, method = "bray")
#make bioenv model against whole meta table with ".", and use gower metric to measure distance in order to incroprate the factor variables. Model will be used up to 7 vriables in final model. parelle processing on 8 cores
#nifH.bioenv <- bioenv(abund_distnifH ~ ., meta3, index = "bray", method = "pearson",
 #                    metric = "gower", upto = 10, parallel = 10)

#summary(nifH.bioenv)
#nifH.bioenvdist<-bioenvdist(nifH.bioenv, which = "best")
#mantel(nifH.bioenvdist, abund_distnifH)
#nifH.bioenv$whichbest


####The best model from the inital BioEnv model selection shows that the Site explains 59% of the variation in the bray-curitis distance. But the best model explains 63% of the variation when Organic Matter and nitritie and nitrate are added in.

Try BioENV without site, this will favor farm managment variables since they are nested within.

Remove site from meta table

meta4<-meta3[,-c(1)]

Rerun BioEnv

#nifH.bioenv.site.na <- bioenv(abund_distnifH ~ ., meta4, index = "bray", 
      #                        method = "pearson", metric = "gower", upto = 10, parallel = 10)

#summary(nifH.bioenv.site.na)
#nifH.bioenvdist.site.na<-bioenvdist(nifH.bioenv.site.na, which = "best")
#mantel(nifH.bioenvdist.site.na, abund_distnifH)
#nifH.bioenv.site.na$whichbest

When Site is removed from the model selection we can see that there is still a bias for the Farm managment variables because they are nested with the site factor. But we do see more come out to build a bigger model. With Tillage + prev_crop + Organic_Matter + Nitrate_Nitrite + Ammonia + Av_Potassium + Manganese + Zinc.


Bioenv with chem data only

#nifH.bioenv.chem <- bioenv(abund_tablenifH ~ Organic_Matter + Moisture_Content + 
 #                               Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + 
  #                              pH + Barium + Calcium + Cobalt + Copper + Iron + Magnesium +
   #                             Manganese + Nickel + Phosphorus + Sulfur + Zinc, meta3, 
    #                          index = "bray", method = "pearson", metric = "gower",
     #                         upto = 10, parallel = 10)

#summary(nifH.bioenv.chem)
#nifH.bioenvdist.chem<-bioenvdist(nifH.bioenv.site.na, which = "best")
#mantel(nifH.bioenvdist.site.na, abund_distnifH)
#nifH.bioenv.chem$whichbest

cca/ordistep model selection

CCA model selection uses a procedure to take the a constrained distance ordination with the complete model and compare it with a unconstrainted model. The model starts with no variables then adds variables that make the best model. This is more

m1_nifH <- cca(abund_tablenifH ~ ., meta3)
m0_nifH <- cca(abund_tablenifH ~ 1, meta3)
m1_nifH
## Call: cca(formula = abund_tablenifH ~ Site + Pea_variety + Plot +
## season_precip + irrgation + total_precip_irr + Tillage + prev_crop
## + grain_yield + Organic_Matter + Moisture_Content +
## Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH +
## Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese
## + Nickel + Phosphorus + Sulfur + Zinc, data = meta3)
## 
##               Inertia Proportion Rank
## Total         10.7659     1.0000     
## Constrained    7.8889     0.7328   33
## Unconstrained  2.8771     0.2672   20
## Inertia is scaled Chi-square 
## Some constraints were aliased because they were collinear (redundant)
## 
## Eigenvalues for constrained axes:
##   CCA1   CCA2   CCA3   CCA4   CCA5   CCA6   CCA7   CCA8   CCA9  CCA10 
## 0.8070 0.6757 0.6117 0.5408 0.5291 0.4814 0.4024 0.3842 0.3332 0.3174 
##  CCA11  CCA12  CCA13  CCA14  CCA15  CCA16  CCA17  CCA18  CCA19  CCA20 
## 0.2820 0.2391 0.2262 0.2006 0.1894 0.1707 0.1638 0.1537 0.1350 0.1302 
##  CCA21  CCA22  CCA23  CCA24  CCA25  CCA26  CCA27  CCA28  CCA29  CCA30 
## 0.1147 0.1077 0.0963 0.0872 0.0804 0.0770 0.0695 0.0613 0.0515 0.0497 
##  CCA31  CCA32  CCA33 
## 0.0441 0.0429 0.0330 
## 
## Eigenvalues for unconstrained axes:
##     CA1     CA2     CA3     CA4     CA5     CA6     CA7     CA8 
## 0.29098 0.28225 0.24918 0.23022 0.20041 0.17909 0.17644 0.16025 
## (Showed only 8 of all 20 unconstrained eigenvalues)
m0_nifH
## Call: cca(formula = abund_tablenifH ~ 1, data = meta3)
## 
##               Inertia Rank
## Total           10.77     
## Unconstrained   10.77   53
## Inertia is scaled Chi-square 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8 
## 0.8189 0.6921 0.6401 0.5860 0.5584 0.5037 0.4818 0.4700 
## (Showed only 8 of all 53 unconstrained eigenvalues)

Ordistep

model_nifH <-ordistep(m0_nifH, scope=formula(m1_nifH))
## 
## Start: abund_tablenifH ~ 1 
## 
##                    Df    AIC      F Pr(>F)   
## + Site              8 868.86 3.0678  0.005 **
## + Tillage           2 873.90 3.2418  0.005 **
## + Organic_Matter    1 874.75 3.5957  0.005 **
## + prev_crop         3 875.14 2.3857  0.005 **
## + Nitrate_Nitrite   1 875.40 2.9283  0.005 **
## + irrgation         1 875.41 2.9174  0.005 **
## + season_precip     1 875.67 2.6597  0.005 **
## + Zinc              1 875.68 2.6463  0.005 **
## + Av_Potassium      1 875.71 2.6178  0.005 **
## + total_precip_irr  1 875.74 2.5881  0.005 **
## + Nickel            1 875.88 2.4512  0.005 **
## + Cobalt            1 875.92 2.4074  0.005 **
## + Barium            1 875.93 2.3929  0.005 **
## + Iron              1 875.93 2.3923  0.005 **
## + Copper            1 875.99 2.3375  0.005 **
## + Plot              1 876.00 2.3303  0.005 **
## + Magnesium         1 876.00 2.3264  0.005 **
## + Av_Phosphorus     1 876.01 2.3202  0.005 **
## + Phosphorus        1 876.07 2.2549  0.005 **
## + Moisture_Content  1 876.18 2.1464  0.005 **
## + Manganese         1 876.18 2.1462  0.005 **
## + grain_yield       1 876.27 2.0533  0.005 **
## + Ammonia           1 876.47 1.8528  0.005 **
## + Sulfur            1 876.47 1.8514  0.005 **
## + Calcium           1 876.57 1.7584  0.005 **
## + pH                1 876.84 1.4897  0.015 * 
## + Pea_variety       6 882.82 0.8462  0.980   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Site 
## 
##        Df    AIC      F Pr(>F)   
## - Site  8 876.36 3.0678  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + Nitrate_Nitrite   1 868.08 2.3207  0.005 **
## + Ammonia           1 869.09 1.4610  0.030 * 
## + Av_Phosphorus     1 869.18 1.3900  0.060 . 
## + Pea_variety       6 872.61 1.0730  0.185   
## + Cobalt            1 869.64 1.0051  0.370   
## + Manganese         1 869.65 0.9978  0.420   
## + grain_yield       1 869.68 0.9722  0.445   
## + Iron              1 869.74 0.9187  0.460   
## + Phosphorus        1 869.79 0.8751  0.565   
## + Copper            1 869.81 0.8624  0.565   
## + Zinc              1 869.77 0.8984  0.575   
## + Barium            1 869.78 0.8907  0.600   
## + Av_Potassium      1 869.86 0.8240  0.705   
## + Magnesium         1 869.87 0.8137  0.705   
## + Nickel            1 869.82 0.8508  0.730   
## + pH                1 869.81 0.8615  0.750   
## + Organic_Matter    1 870.05 0.6648  0.900   
## + Moisture_Content  1 869.98 0.7224  0.910   
## + Calcium           1 870.08 0.6382  0.925   
## + Sulfur            1 870.03 0.6823  0.950   
## + Plot              0 868.86                 
## + season_precip     0 868.86                 
## + irrgation         0 868.86                 
## + total_precip_irr  0 868.86                 
## + Tillage           0 868.86                 
## + prev_crop         0 868.86                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Site + Nitrate_Nitrite 
## 
##                   Df    AIC      F Pr(>F)   
## - Nitrate_Nitrite  1 868.86 2.3207  0.005 **
## - Site             8 875.40 2.9709  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)  
## + Ammonia           1 868.19 1.5305  0.025 *
## + Av_Phosphorus     1 868.30 1.4419  0.070 .
## + Pea_variety       6 871.61 1.0759  0.190  
## + Cobalt            1 868.83 1.0065  0.445  
## + Zinc              1 868.92 0.9313  0.450  
## + Manganese         1 868.81 1.0230  0.455  
## + Iron              1 868.90 0.9530  0.495  
## + Phosphorus        1 868.90 0.9475  0.525  
## + Copper            1 868.95 0.9077  0.550  
## + Barium            1 868.93 0.9309  0.580  
## + Av_Potassium      1 868.97 0.8946  0.600  
## + Magnesium         1 869.03 0.8495  0.665  
## + pH                1 868.99 0.8775  0.690  
## + Nickel            1 869.05 0.8322  0.765  
## + grain_yield       1 869.10 0.7851  0.840  
## + Organic_Matter    1 869.20 0.7088  0.875  
## + Moisture_Content  1 869.15 0.7450  0.900  
## + Sulfur            1 869.21 0.6989  0.920  
## + Calcium           1 869.26 0.6602  0.925  
## + Plot              0 868.08                
## + season_precip     0 868.08                
## + irrgation         0 868.08                
## + total_precip_irr  0 868.08                
## + Tillage           0 868.08                
## + prev_crop         0 868.08                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Site + Nitrate_Nitrite + Ammonia 
## 
##                   Df    AIC      F Pr(>F)   
## - Ammonia          1 868.08 1.5305  0.035 * 
## - Nitrate_Nitrite  1 869.09 2.3726  0.005 **
## - Site             8 875.48 2.8982  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)  
## + Av_Phosphorus     1 868.52 1.3192  0.080 .
## + Pea_variety       6 871.49 1.0780  0.155  
## + Iron              1 868.97 0.9609  0.425  
## + Cobalt            1 868.93 0.9935  0.470  
## + Manganese         1 868.94 0.9858  0.475  
## + Zinc              1 869.01 0.9343  0.505  
## + Barium            1 869.00 0.9394  0.515  
## + Copper            1 869.03 0.9155  0.535  
## + Phosphorus        1 868.98 0.9541  0.545  
## + pH                1 869.03 0.9173  0.620  
## + Magnesium         1 869.10 0.8613  0.670  
## + Av_Potassium      1 869.11 0.8515  0.720  
## + Nickel            1 869.12 0.8419  0.740  
## + grain_yield       1 869.16 0.8129  0.825  
## + Organic_Matter    1 869.28 0.7172  0.885  
## + Sulfur            1 869.27 0.7237  0.905  
## + Calcium           1 869.33 0.6769  0.935  
## + Moisture_Content  1 869.30 0.6981  0.960  
## + Plot              0 868.19                
## + season_precip     0 868.19                
## + irrgation         0 868.19                
## + total_precip_irr  0 868.19                
## + Tillage           0 868.19                
## + prev_crop         0 868.19                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cca without site

m1_nifH_site_na <- cca(abund_tablenifH ~ ., meta4)
m0_nifH_site_na <- cca(abund_tablenifH ~ 1, meta4)
m1_nifH_site_na
## Call: cca(formula = abund_tablenifH ~ Pea_variety + Plot +
## season_precip + irrgation + total_precip_irr + Tillage + prev_crop
## + grain_yield + Organic_Matter + Moisture_Content +
## Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH +
## Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese
## + Nickel + Phosphorus + Sulfur + Zinc, data = meta4)
## 
##               Inertia Proportion Rank
## Total         10.7659     1.0000     
## Constrained    7.6960     0.7148   32
## Unconstrained  3.0700     0.2852   21
## Inertia is scaled Chi-square 
## Some constraints were aliased because they were collinear (redundant)
## 
## Eigenvalues for constrained axes:
##   CCA1   CCA2   CCA3   CCA4   CCA5   CCA6   CCA7   CCA8   CCA9  CCA10 
## 0.8067 0.6725 0.6096 0.5394 0.5055 0.4714 0.3948 0.3835 0.3273 0.3114 
##  CCA11  CCA12  CCA13  CCA14  CCA15  CCA16  CCA17  CCA18  CCA19  CCA20 
## 0.2705 0.2277 0.2256 0.1944 0.1815 0.1687 0.1585 0.1494 0.1316 0.1247 
##  CCA21  CCA22  CCA23  CCA24  CCA25  CCA26  CCA27  CCA28  CCA29  CCA30 
## 0.1138 0.1030 0.0960 0.0871 0.0803 0.0698 0.0633 0.0572 0.0498 0.0445 
##  CCA31  CCA32 
## 0.0430 0.0334 
## 
## Eigenvalues for unconstrained axes:
##     CA1     CA2     CA3     CA4     CA5     CA6     CA7     CA8 
## 0.29949 0.28234 0.26627 0.23089 0.20393 0.19999 0.17647 0.17203 
## (Showed only 8 of all 21 unconstrained eigenvalues)
m0_nifH_site_na
## Call: cca(formula = abund_tablenifH ~ 1, data = meta4)
## 
##               Inertia Rank
## Total           10.77     
## Unconstrained   10.77   53
## Inertia is scaled Chi-square 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8 
## 0.8189 0.6921 0.6401 0.5860 0.5584 0.5037 0.4818 0.4700 
## (Showed only 8 of all 53 unconstrained eigenvalues)

Ordistep

model_nifH <-ordistep(m0_nifH_site_na, scope=formula(m1_nifH_site_na))
## 
## Start: abund_tablenifH ~ 1 
## 
##                    Df    AIC      F Pr(>F)   
## + Tillage           2 873.90 3.2418  0.005 **
## + Organic_Matter    1 874.75 3.5957  0.005 **
## + prev_crop         3 875.14 2.3857  0.005 **
## + Nitrate_Nitrite   1 875.40 2.9283  0.005 **
## + irrgation         1 875.41 2.9174  0.005 **
## + season_precip     1 875.67 2.6597  0.005 **
## + Zinc              1 875.68 2.6463  0.005 **
## + Av_Potassium      1 875.71 2.6178  0.005 **
## + total_precip_irr  1 875.74 2.5881  0.005 **
## + Nickel            1 875.88 2.4512  0.005 **
## + Cobalt            1 875.92 2.4074  0.005 **
## + Barium            1 875.93 2.3929  0.005 **
## + Iron              1 875.93 2.3923  0.005 **
## + Copper            1 875.99 2.3375  0.005 **
## + Plot              1 876.00 2.3303  0.005 **
## + Magnesium         1 876.00 2.3264  0.005 **
## + Av_Phosphorus     1 876.01 2.3202  0.005 **
## + Phosphorus        1 876.07 2.2549  0.005 **
## + Moisture_Content  1 876.18 2.1464  0.005 **
## + Manganese         1 876.18 2.1462  0.005 **
## + grain_yield       1 876.27 2.0533  0.005 **
## + Ammonia           1 876.47 1.8528  0.005 **
## + Calcium           1 876.57 1.7584  0.005 **
## + Sulfur            1 876.47 1.8514  0.010 **
## + pH                1 876.84 1.4897  0.020 * 
## + Pea_variety       6 882.82 0.8462  0.980   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Tillage 
## 
##           Df    AIC      F Pr(>F)   
## - Tillage  2 876.36 3.2418  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + prev_crop         3 871.82 2.5825  0.005 **
## + Organic_Matter    1 872.19 3.5560  0.005 **
## + total_precip_irr  1 872.60 3.1476  0.005 **
## + irrgation         1 872.87 2.8842  0.005 **
## + Av_Potassium      1 872.95 2.8096  0.005 **
## + Av_Phosphorus     1 872.98 2.7771  0.005 **
## + season_precip     1 873.16 2.6012  0.005 **
## + Magnesium         1 873.34 2.4237  0.005 **
## + Phosphorus        1 873.44 2.3340  0.005 **
## + Barium            1 873.62 2.1603  0.005 **
## + Nickel            1 873.63 2.1424  0.005 **
## + Moisture_Content  1 873.65 2.1231  0.005 **
## + Plot              1 873.71 2.0722  0.005 **
## + Manganese         1 873.73 2.0526  0.005 **
## + Nitrate_Nitrite   1 873.77 2.0084  0.005 **
## + Zinc              1 873.94 1.8456  0.005 **
## + Ammonia           1 873.95 1.8361  0.005 **
## + Calcium           1 874.00 1.7954  0.005 **
## + Iron              1 874.06 1.7324  0.005 **
## + Sulfur            1 874.08 1.7112  0.005 **
## + Copper            1 874.13 1.6675  0.005 **
## + grain_yield       1 874.17 1.6277  0.005 **
## + Cobalt            1 873.91 1.8798  0.015 * 
## + pH                1 874.58 1.2350  0.130   
## + Pea_variety       6 879.71 0.9109  0.850   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Tillage + prev_crop 
## 
##             Df    AIC      F Pr(>F)   
## - prev_crop  3 873.90 2.5825  0.005 **
## - Tillage    2 875.14 3.4833  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + irrgation         1 870.32 3.1471  0.005 **
## + Organic_Matter    1 870.67 2.8207  0.005 **
## + Ammonia           1 870.86 2.6513  0.005 **
## + total_precip_irr  1 871.34 2.2116  0.005 **
## + grain_yield       1 871.34 2.2108  0.005 **
## + Magnesium         1 871.34 2.2086  0.005 **
## + Barium            1 871.51 2.0507  0.005 **
## + season_precip     1 871.56 2.0064  0.005 **
## + Nickel            1 871.58 1.9892  0.005 **
## + Manganese         1 871.63 1.9468  0.005 **
## + pH                1 871.84 1.7541  0.005 **
## + Nitrate_Nitrite   1 871.63 1.9486  0.010 **
## + Cobalt            1 871.63 1.9418  0.010 **
## + Av_Potassium      1 871.94 1.6614  0.010 **
## + Calcium           1 872.01 1.6001  0.010 **
## + Zinc              1 871.95 1.6604  0.015 * 
## + Copper            1 872.09 1.5325  0.020 * 
## + Iron              1 871.81 1.7800  0.030 * 
## + Av_Phosphorus     1 872.28 1.3598  0.045 * 
## + Phosphorus        1 872.28 1.3633  0.120   
## + Moisture_Content  1 872.41 1.2435  0.155   
## + Sulfur            1 872.55 1.1177  0.340   
## + Pea_variety       6 876.76 0.9773  0.570   
## + Plot              0 871.82                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation 
## 
##             Df    AIC      F Pr(>F)   
## - irrgation  1 871.82 3.1471  0.005 **
## - prev_crop  3 872.87 2.6883  0.005 **
## - Tillage    2 872.96 3.0756  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + Ammonia           1 870.03 1.9949  0.005 **
## + Manganese         1 870.12 1.9138  0.005 **
## + season_precip     1 870.13 1.9074  0.005 **
## + total_precip_irr  1 870.13 1.9074  0.005 **
## + grain_yield       1 870.20 1.8412  0.005 **
## + pH                1 870.27 1.7801  0.005 **
## + Av_Potassium      1 870.41 1.6562  0.005 **
## + Zinc              1 870.45 1.6182  0.005 **
## + Nitrate_Nitrite   1 870.01 2.0078  0.010 **
## + Organic_Matter    1 870.41 1.6546  0.010 **
## + Nickel            1 870.57 1.5191  0.010 **
## + Magnesium         1 870.54 1.5390  0.030 * 
## + Copper            1 870.59 1.5016  0.040 * 
## + Av_Phosphorus     1 870.70 1.3994  0.060 . 
## + Cobalt            1 870.40 1.6632  0.090 . 
## + Iron              1 870.65 1.4446  0.105   
## + Moisture_Content  1 870.76 1.3526  0.115   
## + Barium            1 870.95 1.1863  0.165   
## + Phosphorus        1 870.83 1.2836  0.220   
## + Sulfur            1 871.08 1.0655  0.375   
## + Pea_variety       6 874.92 1.0036  0.520   
## + Calcium           1 871.27 0.9011  0.685   
## + Plot              0 870.32                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation + Ammonia 
## 
##             Df    AIC      F Pr(>F)   
## - Ammonia    1 870.32 1.9949  0.005 **
## - irrgation  1 870.86 2.4742  0.005 **
## - Tillage    2 872.15 2.7622  0.005 **
## - prev_crop  3 872.66 2.6590  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + season_precip     1 869.61 2.0562  0.005 **
## + total_precip_irr  1 869.61 2.0562  0.005 **
## + Nitrate_Nitrite   1 869.62 2.0545  0.005 **
## + pH                1 869.82 1.8745  0.005 **
## + Manganese         1 869.88 1.8270  0.005 **
## + Av_Potassium      1 870.01 1.7141  0.010 **
## + Organic_Matter    1 870.04 1.6910  0.010 **
## + Zinc              1 870.10 1.6350  0.010 **
## + grain_yield       1 870.11 1.6301  0.010 **
## + Nickel            1 870.13 1.6135  0.020 * 
## + Copper            1 870.20 1.5505  0.025 * 
## + Magnesium         1 870.33 1.4337  0.055 . 
## + Iron              1 870.39 1.3872  0.135   
## + Cobalt            1 870.43 1.3555  0.220   
## + Av_Phosphorus     1 870.75 1.0745  0.250   
## + Pea_variety       6 874.34 1.0205  0.365   
## + Phosphorus        1 870.85 0.9908  0.435   
## + Barium            1 870.85 0.9879  0.510   
## + Sulfur            1 870.95 0.9063  0.625   
## + Calcium           1 870.93 0.9229  0.655   
## + Moisture_Content  1 871.02 0.8480  0.730   
## + Plot              0 870.03                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation + Ammonia +      season_precip 
## 
##                 Df    AIC      F Pr(>F)   
## - season_precip  1 870.03 2.0562  0.005 **
## - Ammonia        1 870.13 2.1422  0.005 **
## - irrgation      1 870.58 2.5391  0.005 **
## - Tillage        2 871.49 2.5867  0.005 **
## - prev_crop      3 871.65 2.4051  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + grain_yield       1 869.54 1.7222  0.010 **
## + Nitrate_Nitrite   1 868.86 2.3031  0.020 * 
## + Cobalt            1 869.93 1.3919  0.160   
## + Manganese         1 869.86 1.4507  0.170   
## + Av_Phosphorus     1 870.26 1.1212  0.270   
## + Pea_variety       6 873.52 1.0515  0.310   
## + Magnesium         1 870.30 1.0830  0.355   
## + Iron              1 870.31 1.0795  0.360   
## + pH                1 870.37 1.0279  0.415   
## + Av_Potassium      1 870.39 1.0079  0.460   
## + Phosphorus        1 870.39 1.0080  0.470   
## + Barium            1 870.41 0.9929  0.515   
## + Zinc              1 870.53 0.8910  0.575   
## + Organic_Matter    1 870.50 0.9177  0.610   
## + Nickel            1 870.60 0.8379  0.645   
## + Sulfur            1 870.54 0.8863  0.660   
## + Copper            1 870.61 0.8285  0.680   
## + Moisture_Content  1 870.68 0.7700  0.800   
## + Calcium           1 870.82 0.6544  0.915   
## + Plot              0 869.61                 
## + total_precip_irr  0 869.61                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation + Ammonia +      season_precip + grain_yield 
## 
##                 Df    AIC      F Pr(>F)   
## - Ammonia        1 869.37 1.5161  0.015 * 
## - grain_yield    1 869.61 1.7222  0.005 **
## - season_precip  1 870.11 2.1400  0.005 **
## - irrgation      1 870.60 2.5682  0.005 **
## - Tillage        2 871.60 2.6131  0.005 **
## - prev_crop      3 871.71 2.3952  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + Nitrate_Nitrite   1 868.62 2.3907  0.010 **
## + Av_Phosphorus     1 870.05 1.2029  0.160   
## + Pea_variety       6 873.33 1.0408  0.290   
## + Cobalt            1 870.30 0.9985  0.415   
## + Av_Potassium      1 870.34 0.9699  0.420   
## + Iron              1 870.36 0.9536  0.425   
## + Manganese         1 870.25 1.0415  0.465   
## + pH                1 870.30 1.0004  0.505   
## + Phosphorus        1 870.35 0.9605  0.535   
## + Zinc              1 870.44 0.8877  0.540   
## + Magnesium         1 870.46 0.8682  0.605   
## + Barium            1 870.42 0.9015  0.620   
## + Copper            1 870.49 0.8450  0.685   
## + Nickel            1 870.49 0.8494  0.740   
## + Sulfur            1 870.62 0.7421  0.830   
## + Organic_Matter    1 870.62 0.7391  0.850   
## + Moisture_Content  1 870.65 0.7188  0.890   
## + Calcium           1 870.70 0.6781  0.905   
## + Plot              0 869.54                 
## + total_precip_irr  0 869.54                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Tillage + prev_crop + irrgation + Ammonia +      season_precip + grain_yield + Nitrate_Nitrite 
## 
##                   Df    AIC      F Pr(>F)   
## - Ammonia          1 868.61 1.6173  0.015 * 
## - grain_yield      1 868.86 1.8213  0.005 **
## - season_precip    1 869.33 2.2157  0.005 **
## - Nitrate_Nitrite  1 869.54 2.3907  0.005 **
## - irrgation        1 869.85 2.6529  0.005 **
## - Tillage          2 870.93 2.6640  0.005 **
## - prev_crop        3 871.16 2.4551  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)
## + Av_Phosphorus     1 869.14 1.1639  0.225
## + Pea_variety       6 872.23 1.0371  0.365
## + Cobalt            1 869.31 1.0309  0.370
## + Manganese         1 869.27 1.0642  0.380
## + Av_Potassium      1 869.31 1.0299  0.390
## + pH                1 869.33 1.0128  0.435
## + Phosphorus        1 869.31 1.0293  0.450
## + Iron              1 869.37 0.9871  0.485
## + Zinc              1 869.47 0.9046  0.525
## + Magnesium         1 869.45 0.9228  0.580
## + Barium            1 869.46 0.9115  0.595
## + Copper            1 869.46 0.9110  0.605
## + Nickel            1 869.55 0.8407  0.780
## + Moisture_Content  1 869.71 0.7150  0.880
## + Calcium           1 869.73 0.6989  0.880
## + Sulfur            1 869.67 0.7431  0.890
## + Organic_Matter    1 869.68 0.7370  0.905
## + Plot              0 868.62              
## + total_precip_irr  0 868.62
m1_nifH_cca_chem <- cca(abund_tablenifH ~ Organic_Matter + Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium +  pH + Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc , meta3)
m0_nifH_cca_chem <- cca(abund_tablenifH ~ 1, meta3)
m1_nifH_cca_chem
## Call: cca(formula = abund_tablenifH ~ Organic_Matter +
## Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus +
## Av_Potassium + pH + Barium + Calcium + Cobalt + Copper + Iron +
## Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc, data
## = meta3)
## 
##               Inertia Proportion Rank
## Total         10.7659     1.0000     
## Constrained    5.2066     0.4836   18
## Unconstrained  5.5593     0.5164   35
## Inertia is scaled Chi-square 
## 
## Eigenvalues for constrained axes:
##   CCA1   CCA2   CCA3   CCA4   CCA5   CCA6   CCA7   CCA8   CCA9  CCA10 
## 0.7916 0.6368 0.5430 0.5237 0.4736 0.4137 0.3403 0.2283 0.2008 0.1907 
##  CCA11  CCA12  CCA13  CCA14  CCA15  CCA16  CCA17  CCA18 
## 0.1609 0.1430 0.1267 0.1085 0.0994 0.0838 0.0794 0.0625 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8 
## 0.4573 0.4358 0.3745 0.3252 0.2860 0.2803 0.2665 0.2449 
## (Showed only 8 of all 35 unconstrained eigenvalues)
m0_nifH_cca_chem
## Call: cca(formula = abund_tablenifH ~ 1, data = meta3)
## 
##               Inertia Rank
## Total           10.77     
## Unconstrained   10.77   53
## Inertia is scaled Chi-square 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8 
## 0.8189 0.6921 0.6401 0.5860 0.5584 0.5037 0.4818 0.4700 
## (Showed only 8 of all 53 unconstrained eigenvalues)
model_nifH_cca_chem <-ordistep(m0_nifH_cca_chem, scope=formula(m1_nifH_cca_chem))
## 
## Start: abund_tablenifH ~ 1 
## 
##                    Df    AIC      F Pr(>F)   
## + Organic_Matter    1 874.75 3.5957  0.005 **
## + Nitrate_Nitrite   1 875.40 2.9283  0.005 **
## + Zinc              1 875.68 2.6463  0.005 **
## + Av_Potassium      1 875.71 2.6178  0.005 **
## + Nickel            1 875.88 2.4512  0.005 **
## + Cobalt            1 875.92 2.4074  0.005 **
## + Barium            1 875.93 2.3929  0.005 **
## + Iron              1 875.93 2.3923  0.005 **
## + Copper            1 875.99 2.3375  0.005 **
## + Magnesium         1 876.00 2.3264  0.005 **
## + Av_Phosphorus     1 876.01 2.3202  0.005 **
## + Phosphorus        1 876.07 2.2549  0.005 **
## + Moisture_Content  1 876.18 2.1464  0.005 **
## + Manganese         1 876.18 2.1462  0.005 **
## + Ammonia           1 876.47 1.8528  0.005 **
## + Calcium           1 876.57 1.7584  0.005 **
## + Sulfur            1 876.47 1.8514  0.010 **
## + pH                1 876.84 1.4897  0.030 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Organic_Matter 
## 
##                  Df    AIC      F Pr(>F)   
## - Organic_Matter  1 876.36 3.5957  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + Nitrate_Nitrite   1 873.56 3.1042  0.005 **
## + Zinc              1 873.97 2.6966  0.005 **
## + Barium            1 874.16 2.5125  0.005 **
## + Magnesium         1 874.20 2.4698  0.005 **
## + Cobalt            1 874.21 2.4614  0.005 **
## + Copper            1 874.22 2.4475  0.005 **
## + Av_Potassium      1 874.23 2.4431  0.005 **
## + Iron              1 874.24 2.4244  0.005 **
## + Av_Phosphorus     1 874.27 2.3965  0.005 **
## + Phosphorus        1 874.34 2.3325  0.005 **
## + Nickel            1 874.36 2.3098  0.005 **
## + Manganese         1 874.39 2.2789  0.005 **
## + Moisture_Content  1 874.48 2.1872  0.005 **
## + Sulfur            1 874.63 2.0475  0.005 **
## + Calcium           1 874.79 1.8862  0.005 **
## + Ammonia           1 874.88 1.7998  0.005 **
## + pH                1 875.14 1.5430  0.010 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite 
## 
##                   Df    AIC      F Pr(>F)   
## - Nitrate_Nitrite  1 874.75 3.1042  0.005 **
## - Organic_Matter   1 875.40 3.7617  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + Av_Potassium      1 872.79 2.6309  0.005 **
## + Magnesium         1 872.92 2.5090  0.005 **
## + Av_Phosphorus     1 872.98 2.4507  0.005 **
## + Phosphorus        1 872.99 2.4422  0.005 **
## + Moisture_Content  1 873.15 2.2870  0.005 **
## + Manganese         1 873.22 2.2192  0.005 **
## + Zinc              1 873.36 2.0786  0.005 **
## + Sulfur            1 873.39 2.0553  0.005 **
## + Copper            1 873.50 1.9447  0.005 **
## + Cobalt            1 873.54 1.9069  0.005 **
## + Calcium           1 873.55 1.8973  0.005 **
## + Barium            1 873.61 1.8369  0.005 **
## + Ammonia           1 873.65 1.8007  0.005 **
## + Iron              1 873.67 1.7782  0.005 **
## + Nickel            1 873.84 1.6164  0.005 **
## + pH                1 873.98 1.4887  0.020 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium 
## 
##                   Df    AIC      F Pr(>F)   
## - Av_Potassium     1 873.56 2.6309  0.005 **
## - Nitrate_Nitrite  1 874.23 3.2820  0.005 **
## - Organic_Matter   1 874.51 3.5682  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + Av_Phosphorus     1 872.25 2.3577  0.005 **
## + Magnesium         1 872.36 2.2606  0.005 **
## + Manganese         1 872.56 2.0639  0.005 **
## + Sulfur            1 872.61 2.0249  0.005 **
## + Calcium           1 872.63 2.0019  0.005 **
## + Ammonia           1 872.75 1.8850  0.005 **
## + Barium            1 872.80 1.8420  0.005 **
## + Nickel            1 872.90 1.7454  0.005 **
## + pH                1 872.93 1.7153  0.005 **
## + Copper            1 872.96 1.6867  0.005 **
## + Zinc              1 872.98 1.6722  0.005 **
## + Phosphorus        1 873.06 1.5945  0.005 **
## + Cobalt            1 872.89 1.7581  0.010 **
## + Moisture_Content  1 873.18 1.4875  0.015 * 
## + Iron              1 873.16 1.4996  0.035 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium +      Av_Phosphorus 
## 
##                   Df    AIC      F Pr(>F)   
## - Av_Phosphorus    1 872.79 2.3577  0.005 **
## - Av_Potassium     1 872.98 2.5341  0.005 **
## - Nitrate_Nitrite  1 873.69 3.2208  0.005 **
## - Organic_Matter   1 874.15 3.6662  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + Manganese         1 871.90 2.1402  0.005 **
## + Magnesium         1 872.03 2.0192  0.005 **
## + Sulfur            1 872.11 1.9489  0.005 **
## + Calcium           1 872.19 1.8727  0.005 **
## + Barium            1 872.29 1.7817  0.005 **
## + Ammonia           1 872.36 1.7145  0.005 **
## + Zinc              1 872.38 1.6967  0.005 **
## + Copper            1 872.44 1.6399  0.005 **
## + Phosphorus        1 872.48 1.6044  0.005 **
## + Cobalt            1 872.28 1.7883  0.010 **
## + Nickel            1 872.47 1.6161  0.010 **
## + Iron              1 872.54 1.5510  0.010 **
## + pH                1 872.58 1.5087  0.025 * 
## + Moisture_Content  1 872.82 1.2883  0.080 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium +      Av_Phosphorus + Manganese 
## 
##                   Df    AIC      F Pr(>F)   
## - Manganese        1 872.25 2.1402  0.005 **
## - Av_Phosphorus    1 872.56 2.4287  0.005 **
## - Av_Potassium     1 872.74 2.5890  0.005 **
## - Nitrate_Nitrite  1 873.12 2.9469  0.005 **
## - Organic_Matter   1 873.85 3.6433  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + Barium            1 871.44 2.1898  0.005 **
## + Sulfur            1 871.71 1.9435  0.005 **
## + Magnesium         1 871.76 1.9029  0.005 **
## + Zinc              1 871.79 1.8710  0.005 **
## + Cobalt            1 871.82 1.8413  0.005 **
## + Iron              1 871.91 1.7662  0.005 **
## + Copper            1 871.91 1.7627  0.005 **
## + Calcium           1 871.92 1.7571  0.005 **
## + Nickel            1 872.03 1.6558  0.005 **
## + Phosphorus        1 872.05 1.6338  0.005 **
## + pH                1 872.11 1.5817  0.005 **
## + Ammonia           1 872.20 1.5044  0.010 **
## + Moisture_Content  1 872.38 1.3407  0.050 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium +      Av_Phosphorus + Manganese + Barium 
## 
##                   Df    AIC      F Pr(>F)   
## - Barium           1 871.90 2.1898  0.005 **
## - Av_Phosphorus    1 872.08 2.3589  0.005 **
## - Nitrate_Nitrite  1 872.21 2.4735  0.005 **
## - Av_Potassium     1 872.29 2.5438  0.005 **
## - Manganese        1 872.29 2.5440  0.005 **
## - Organic_Matter   1 873.10 3.2928  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + Cobalt            1 871.21 1.9353  0.005 **
## + Copper            1 871.30 1.8613  0.005 **
## + Sulfur            1 871.30 1.8589  0.005 **
## + Iron              1 871.31 1.8503  0.005 **
## + Magnesium         1 871.35 1.8125  0.005 **
## + Zinc              1 871.42 1.7495  0.005 **
## + Calcium           1 871.43 1.7405  0.005 **
## + Phosphorus        1 871.56 1.6326  0.005 **
## + Nickel            1 871.65 1.5524  0.005 **
## + pH                1 871.56 1.6330  0.010 **
## + Moisture_Content  1 871.84 1.3860  0.015 * 
## + Ammonia           1 871.84 1.3831  0.025 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium +      Av_Phosphorus + Manganese + Barium + Cobalt 
## 
##                   Df    AIC      F Pr(>F)   
## - Cobalt           1 871.44 1.9353  0.005 **
## - Barium           1 871.82 2.2773  0.005 **
## - Av_Phosphorus    1 871.87 2.3225  0.005 **
## - Nitrate_Nitrite  1 871.91 2.3545  0.005 **
## - Manganese        1 871.92 2.3678  0.005 **
## - Av_Potassium     1 872.23 2.6409  0.005 **
## - Organic_Matter   1 873.08 3.4110  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + Sulfur            1 870.86 2.0053  0.005 **
## + Calcium           1 870.96 1.9217  0.005 **
## + Nickel            1 871.19 1.7216  0.005 **
## + pH                1 871.48 1.4699  0.010 **
## + Phosphorus        1 871.57 1.3904  0.025 * 
## + Moisture_Content  1 871.58 1.3814  0.025 * 
## + Ammonia           1 871.61 1.3613  0.040 * 
## + Zinc              1 871.68 1.2995  0.060 . 
## + Magnesium         1 871.71 1.2731  0.070 . 
## + Copper            1 871.87 1.1341  0.215   
## + Iron              1 872.01 1.0150  0.530   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium +      Av_Phosphorus + Manganese + Barium + Cobalt + Sulfur 
## 
##                   Df    AIC      F Pr(>F)   
## - Sulfur           1 871.21 2.0053  0.005 **
## - Cobalt           1 871.30 2.0804  0.005 **
## - Barium           1 871.39 2.1625  0.005 **
## - Av_Phosphorus    1 871.44 2.1990  0.005 **
## - Manganese        1 871.57 2.3147  0.005 **
## - Nitrate_Nitrite  1 871.67 2.4073  0.005 **
## - Av_Potassium     1 872.03 2.7182  0.005 **
## - Organic_Matter   1 872.99 3.5736  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)   
## + Magnesium         1 871.12 1.4435  0.010 **
## + Nickel            1 871.12 1.4399  0.040 * 
## + pH                1 871.36 1.2417  0.090 . 
## + Copper            1 871.38 1.2255  0.125   
## + Phosphorus        1 871.35 1.2508  0.130   
## + Zinc              1 871.41 1.2013  0.145   
## + Ammonia           1 871.53 1.0985  0.305   
## + Moisture_Content  1 871.63 1.0163  0.485   
## + Iron              1 871.66 0.9862  0.520   
## + Calcium           1 871.93 0.7671  0.835   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium +      Av_Phosphorus + Manganese + Barium + Cobalt + Sulfur + Magnesium 
## 
##                   Df    AIC      F Pr(>F)   
## - Cobalt           1 870.55 1.1831  0.225   
## - Magnesium        1 870.86 1.4435  0.030 * 
## - Av_Phosphorus    1 870.94 1.5081  0.005 **
## - Barium           1 871.20 1.7285  0.005 **
## - Sulfur           1 871.71 2.1626  0.005 **
## - Manganese        1 871.72 2.1755  0.005 **
## - Nitrate_Nitrite  1 871.80 2.2421  0.005 **
## - Organic_Matter   1 871.94 2.3606  0.005 **
## - Av_Potassium     1 872.38 2.7380  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium +      Av_Phosphorus + Manganese + Barium + Sulfur + Magnesium 
## 
##                    Df    AIC      F Pr(>F)   
## + Phosphorus        1 870.39 1.7956  0.005 **
## + Iron              1 870.99 1.2859  0.155   
## + Ammonia           1 871.17 1.1384  0.210   
## + pH                1 871.19 1.1215  0.230   
## + Cobalt            1 871.12 1.1831  0.245   
## + Nickel            1 871.20 1.1160  0.310   
## + Zinc              1 871.44 0.9148  0.650   
## + Moisture_Content  1 871.48 0.8821  0.675   
## + Copper            1 871.45 0.9088  0.710   
## + Calcium           1 871.50 0.8643  0.740   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: abund_tablenifH ~ Organic_Matter + Nitrate_Nitrite + Av_Potassium +      Av_Phosphorus + Manganese + Barium + Sulfur + Magnesium +      Phosphorus 
## 
##                   Df    AIC      F Pr(>F)   
## - Av_Phosphorus    1 870.15 1.4610  0.025 * 
## - Av_Potassium     1 870.27 1.5617  0.010 **
## - Phosphorus       1 870.55 1.7956  0.005 **
## - Magnesium        1 870.72 1.9382  0.005 **
## - Nitrate_Nitrite  1 870.93 2.1213  0.005 **
## - Sulfur           1 871.03 2.2060  0.005 **
## - Manganese        1 871.07 2.2387  0.005 **
## - Barium           1 871.10 2.2666  0.005 **
## - Organic_Matter   1 871.48 2.5897  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    Df    AIC      F Pr(>F)
## + pH                1 870.92 1.1871  0.180
## + Ammonia           1 870.94 1.1697  0.215
## + Nickel            1 871.03 1.0927  0.275
## + Cobalt            1 871.04 1.0916  0.275
## + Iron              1 871.03 1.0935  0.315
## + Zinc              1 871.13 1.0140  0.445
## + Calcium           1 871.23 0.9362  0.510
## + Moisture_Content  1 871.28 0.8899  0.735
## + Copper            1 871.33 0.8507  0.775
model_nifH_cca_chem$anova
##                   Df    AIC      F Pr(>F)   
## + Organic_Matter   1 874.75 3.5957  0.005 **
## + Nitrate_Nitrite  1 873.56 3.1042  0.005 **
## + Av_Potassium     1 872.79 2.6309  0.005 **
## + Av_Phosphorus    1 872.25 2.3577  0.005 **
## + Manganese        1 871.90 2.1402  0.005 **
## + Barium           1 871.44 2.1898  0.005 **
## + Cobalt           1 871.21 1.9353  0.005 **
## + Sulfur           1 870.86 2.0053  0.005 **
## + Magnesium        1 871.12 1.4435  0.010 **
## - Cobalt           1 870.55 1.1831  0.225   
## + Phosphorus       1 870.39 1.7956  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

RDA is a linear cca

m1_nifH_rda <- rda(abund_tablenifH ~ ., meta3)
m0_nifH_rda <- rda(abund_tablenifH ~ 1, meta3)
m1_nifH_rda
## Call: rda(formula = abund_tablenifH ~ Site + Pea_variety + Plot +
## season_precip + irrgation + total_precip_irr + Tillage + prev_crop
## + grain_yield + Organic_Matter + Moisture_Content +
## Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium + pH +
## Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese
## + Nickel + Phosphorus + Sulfur + Zinc, data = meta3)
## 
##                 Inertia Proportion Rank
## Total         1.258e+11  1.000e+00     
## Constrained   9.747e+10  7.745e-01   33
## Unconstrained 2.837e+10  2.255e-01   20
## Inertia is variance 
## Some constraints were aliased because they were collinear (redundant)
## 
## Eigenvalues for constrained axes:
##        RDA1        RDA2        RDA3        RDA4        RDA5        RDA6 
## 16353665417 15356453441 12209872418 10446884520  7322044854  5776558730 
##        RDA7        RDA8        RDA9       RDA10       RDA11       RDA12 
##  4833284412  4476474731  3654563170  3048908811  2374942944  1890413915 
##       RDA13       RDA14       RDA15       RDA16       RDA17       RDA18 
##  1488069677  1072841995   989712375   941967373   836809708   626843871 
##       RDA19       RDA20       RDA21       RDA22       RDA23       RDA24 
##   573856607   508536865   403729752   376714333   325054973   282465322 
##       RDA25       RDA26       RDA27       RDA28       RDA29       RDA30 
##   238159099   207287067   185699351   163253845   127967489   116333192 
##       RDA31       RDA32       RDA33 
##    95605085    84351951    75938358 
## 
## Eigenvalues for unconstrained axes:
##        PC1        PC2        PC3        PC4        PC5        PC6 
## 5599447125 4342938458 3880332013 3368259818 2027561149 1483067201 
##        PC7        PC8 
## 1310478544 1156592863 
## (Showed only 8 of all 20 unconstrained eigenvalues)
m0_nifH_rda
## Call: rda(formula = abund_tablenifH ~ 1, data = meta3)
## 
##                 Inertia Rank
## Total         1.258e+11     
## Unconstrained 1.258e+11   53
## Inertia is variance 
## 
## Eigenvalues for unconstrained axes:
##         PC1         PC2         PC3         PC4         PC5         PC6 
## 18388165460 17017988044 13929669751 13228540294 10298920752  7474983345 
##         PC7         PC8 
##  6353906098  5386444205 
## (Showed only 8 of all 53 unconstrained eigenvalues)
plot(m1_nifH_rda)

model_rda_nifH <-ordiR2step(m0_nifH_rda, scope=formula(m1_nifH_rda))
## Step: R2.adj= 0 
## Call: abund_tablenifH ~ 1 
##  
##                    R2.adjusted
## <All variables>     0.40253072
## + Site              0.39284983
## + prev_crop         0.14946826
## + Tillage           0.09118753
## + total_precip_irr  0.06245835
## + Av_Phosphorus     0.06052058
## + Phosphorus        0.05195778
## + Av_Potassium      0.05112278
## + Magnesium         0.04894792
## + irrgation         0.04678398
## + Plot              0.04077308
## + Moisture_Content  0.04047656
## + season_precip     0.03754349
## + Barium            0.03466802
## + Manganese         0.03434338
## + Zinc              0.03212237
## + Organic_Matter    0.03177594
## + Ammonia           0.02812448
## + Nitrate_Nitrite   0.02759831
## + Nickel            0.02460517
## + Copper            0.02201680
## + Sulfur            0.02175055
## + Calcium           0.01980531
## + grain_yield       0.01948985
## + Iron              0.01591920
## + pH                0.01585647
## + Cobalt            0.01268035
## <none>              0.00000000
## + Pea_variety      -0.02472576
## 
##        Df    AIC      F Pr(>F)   
## + Site  8 1361.4 5.2866  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.3928498 
## Call: abund_tablenifH ~ Site 
##  
##                    R2.adjusted
## + Pea_variety        0.4053417
## <All variables>      0.4025307
## + Ammonia            0.3971916
## + Nitrate_Nitrite    0.3945793
## <none>               0.3928498
## + Plot               0.3928498
## + season_precip      0.3928498
## + irrgation          0.3928498
## + total_precip_irr   0.3928498
## + Tillage            0.3928498
## + prev_crop          0.3928498
## + Av_Potassium       0.3920019
## + Organic_Matter     0.3905573
## + Av_Phosphorus      0.3888507
## + Sulfur             0.3888353
## + pH                 0.3883561
## + Calcium            0.3881788
## + Moisture_Content   0.3878418
## + Manganese          0.3866273
## + grain_yield        0.3860617
## + Magnesium          0.3853474
## + Nickel             0.3852829
## + Cobalt             0.3846202
## + Copper             0.3845997
## + Iron               0.3842913
## + Barium             0.3842156
## + Phosphorus         0.3841708
## + Zinc               0.3838492

Chemistry only model

m1_nifH_rda_chem <- rda(abund_tablenifH ~ Organic_Matter + Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium +  pH + Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc , meta3)
m0_nifH_rda_chem <- rda(abund_tablenifH ~ 1, meta3)
m1_nifH_rda_chem
## Call: rda(formula = abund_tablenifH ~ Organic_Matter +
## Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus +
## Av_Potassium + pH + Barium + Calcium + Cobalt + Copper + Iron +
## Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc, data
## = meta3)
## 
##                 Inertia Proportion Rank
## Total         1.258e+11  1.000e+00     
## Constrained   6.893e+10  5.478e-01   18
## Unconstrained 5.690e+10  4.522e-01   35
## Inertia is variance 
## 
## Eigenvalues for constrained axes:
##        RDA1        RDA2        RDA3        RDA4        RDA5        RDA6 
## 14804112412 13108358521 10437172065  8136887370  5686929136  3911450523 
##        RDA7        RDA8        RDA9       RDA10       RDA11       RDA12 
##  3487763724  2618476736  1896564759  1302627898   875978909   653782995 
##       RDA13       RDA14       RDA15       RDA16       RDA17       RDA18 
##   566194235   425743412   312506701   290536181   223053545   193225475 
## 
## Eigenvalues for unconstrained axes:
##         PC1         PC2         PC3         PC4         PC5         PC6 
## 11018590525  7155476437  5939185648  5422891779  4017288576  3016129382 
##         PC7         PC8 
##  2886770730  2071617511 
## (Showed only 8 of all 35 unconstrained eigenvalues)
m0_nifH_rda_chem
## Call: rda(formula = abund_tablenifH ~ 1, data = meta3)
## 
##                 Inertia Rank
## Total         1.258e+11     
## Unconstrained 1.258e+11   53
## Inertia is variance 
## 
## Eigenvalues for unconstrained axes:
##         PC1         PC2         PC3         PC4         PC5         PC6 
## 18388165460 17017988044 13929669751 13228540294 10298920752  7474983345 
##         PC7         PC8 
##  6353906098  5386444205 
## (Showed only 8 of all 53 unconstrained eigenvalues)
plot(m1_nifH_rda_chem)

model_rda_chem_nifH <-ordiR2step(m0_nifH_rda_chem, scope=formula(m1_nifH_rda_chem))
## Step: R2.adj= 0 
## Call: abund_tablenifH ~ 1 
##  
##                    R2.adjusted
## <All variables>     0.31521855
## + Av_Phosphorus     0.06052058
## + Phosphorus        0.05195778
## + Av_Potassium      0.05112278
## + Magnesium         0.04894792
## + Moisture_Content  0.04047656
## + Barium            0.03466802
## + Manganese         0.03434338
## + Zinc              0.03212237
## + Organic_Matter    0.03177594
## + Ammonia           0.02812448
## + Nitrate_Nitrite   0.02759831
## + Nickel            0.02460517
## + Copper            0.02201680
## + Sulfur            0.02175055
## + Calcium           0.01980531
## + Iron              0.01591920
## + pH                0.01585647
## + Cobalt            0.01268035
## <none>              0.00000000
## 
##                 Df    AIC      F Pr(>F)   
## + Av_Phosphorus  1 1378.7 4.4142  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.06052058 
## Call: abund_tablenifH ~ Av_Phosphorus 
##  
##                    R2.adjusted
## <All variables>     0.31521855
## + Av_Potassium      0.10253413
## + Phosphorus        0.09795318
## + Zinc              0.09505177
## + Manganese         0.09455271
## + Organic_Matter    0.09324882
## + Ammonia           0.09270795
## + Nickel            0.08781964
## + Magnesium         0.08490752
## + Barium            0.08424954
## + Moisture_Content  0.08392402
## + Nitrate_Nitrite   0.08340692
## + Sulfur            0.08244671
## + Copper            0.08171114
## + Iron              0.08134433
## + Calcium           0.07884137
## + Cobalt            0.07764768
## + pH                0.07328876
## <none>              0.06052058
## 
##                Df    AIC      F Pr(>F)   
## + Av_Potassium  1 1377.2 3.4343  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.1025341 
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium 
##  
##                    R2.adjusted
## <All variables>      0.3152185
## + Organic_Matter     0.1426010
## + Zinc               0.1346765
## + Manganese          0.1346442
## + Sulfur             0.1321684
## + Nickel             0.1314652
## + Phosphorus         0.1297865
## + Calcium            0.1276961
## + Magnesium          0.1269648
## + Nitrate_Nitrite    0.1268511
## + Ammonia            0.1264398
## + Barium             0.1263189
## + Copper             0.1232790
## + Moisture_Content   0.1227996
## + pH                 0.1220798
## + Iron               0.1192709
## + Cobalt             0.1140726
## <none>               0.1025341
## 
##                  Df    AIC      F Pr(>F)   
## + Organic_Matter  1 1375.7 3.3833  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.142601 
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter 
##  
##                    R2.adjusted
## <All variables>      0.3152185
## + Barium             0.1806009
## + Manganese          0.1778531
## + Zinc               0.1739765
## + Ammonia            0.1728053
## + Sulfur             0.1698627
## + Calcium            0.1691275
## + Magnesium          0.1685205
## + Nitrate_Nitrite    0.1637934
## + Nickel             0.1637017
## + Copper             0.1636724
## + pH                 0.1635979
## + Phosphorus         0.1610160
## + Iron               0.1572052
## + Moisture_Content   0.1536428
## + Cobalt             0.1526631
## <none>               0.1426010
## 
##          Df    AIC      F Pr(>F)   
## + Barium  1 1374.1 3.3188  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.1806009 
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter +      Barium 
##  
##                    R2.adjusted
## <All variables>      0.3152185
## + Manganese          0.2425367
## + Sulfur             0.2152786
## + Calcium            0.2146288
## + Magnesium          0.2075492
## + Cobalt             0.2059481
## + Iron               0.2045766
## + Zinc               0.2041698
## + Nickel             0.2040281
## + pH                 0.2023023
## + Ammonia            0.2016233
## + Copper             0.2015660
## + Phosphorus         0.2011897
## + Moisture_Content   0.1920917
## + Nitrate_Nitrite    0.1880401
## <none>               0.1806009
## 
##             Df    AIC      F Pr(>F)   
## + Manganese  1 1370.8 5.0066  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.2425367 
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter +      Barium + Manganese 
##  
##                    R2.adjusted
## <All variables>      0.3152185
## + Iron               0.2741671
## + Magnesium          0.2715470
## + Zinc               0.2703388
## + Sulfur             0.2688961
## + Copper             0.2683106
## + Cobalt             0.2675120
## + pH                 0.2661270
## + Calcium            0.2650136
## + Phosphorus         0.2645607
## + Nickel             0.2626842
## + Moisture_Content   0.2555896
## + Ammonia            0.2539120
## + Nitrate_Nitrite    0.2511329
## <none>               0.2425367
## 
##        Df    AIC      F Pr(>F)   
## + Iron  1 1369.3 3.0917  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.2741671 
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter +      Barium + Manganese + Iron 
##  
##                    R2.adjusted
## <All variables>      0.3152185
## + Sulfur             0.3077588
## + Calcium            0.3034844
## + pH                 0.2895106
## + Phosphorus         0.2873660
## + Moisture_Content   0.2869075
## + Ammonia            0.2867278
## + Nickel             0.2850464
## + Zinc               0.2833410
## + Magnesium          0.2802296
## + Nitrate_Nitrite    0.2787387
## + Cobalt             0.2743731
## <none>               0.2741671
## + Copper             0.2741048
## 
##          Df    AIC      F Pr(>F)   
## + Sulfur  1 1367.6 3.2807  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.3077588 
## Call: abund_tablenifH ~ Av_Phosphorus + Av_Potassium + Organic_Matter +      Barium + Manganese + Iron + Sulfur 
##  
##                    R2.adjusted
## + Magnesium          0.3193868
## + Zinc               0.3191464
## <All variables>      0.3152185
## + Ammonia            0.3141646
## + pH                 0.3132794
## + Phosphorus         0.3131892
## + Nitrate_Nitrite    0.3129262
## + Copper             0.3093376
## + Moisture_Content   0.3090885
## + Nickel             0.3081999
## <none>               0.3077588
## + Cobalt             0.3049041
## + Calcium            0.3007448
model_rda_chem_nifH$anova
##                    R2.adj Df    AIC      F Pr(>F)   
## + Av_Phosphorus  0.060521  1 1378.7 4.4142  0.002 **
## + Av_Potassium   0.102534  1 1377.2 3.4343  0.002 **
## + Organic_Matter 0.142601  1 1375.7 3.3833  0.002 **
## + Barium         0.180601  1 1374.1 3.3188  0.002 **
## + Manganese      0.242537  1 1370.8 5.0066  0.002 **
## + Iron           0.274167  1 1369.3 3.0917  0.002 **
## + Sulfur         0.307759  1 1367.6 3.2807  0.002 **
## <All variables>  0.315219                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Our RDA chem model can explain 30% of the variance in the bacterial community.

CAP Ordination model building

m1_nifH_cap_chem<- capscale(abund_tablenifH ~ Organic_Matter + Moisture_Content + Nitrate_Nitrite + Ammonia + Av_Phosphorus + Av_Potassium +  pH + Barium + Calcium + Cobalt + Copper + Iron + Magnesium + Manganese + Nickel + Phosphorus + Sulfur + Zinc , data = meta3, distance = "bray")
m0_nifH_cap_chem<- capscale(abund_distnifH ~ 1, meta3)
plot(m1_nifH_cap_chem)

model_cap_chem_nifH <-ordiR2step(m0_nifH_cap_chem, scope=formula(m1_nifH_cap_chem))
## Step: R2.adj= 0 
## Call: abund_distnifH ~ 1 
##  
##                    R2.adjusted
## <All variables>     0.42412042
## + Organic_Matter    0.08664823
## + Av_Potassium      0.06424058
## + Av_Phosphorus     0.05197708
## + Zinc              0.04843832
## + Nitrate_Nitrite   0.04757351
## + Nickel            0.04703479
## + Barium            0.04533845
## + Iron              0.04135991
## + Magnesium         0.04111715
## + Cobalt            0.03937854
## + Phosphorus        0.03739460
## + Manganese         0.03679948
## + Copper            0.03582041
## + Ammonia           0.03314279
## + Moisture_Content  0.02637968
## + Sulfur            0.02145207
## + pH                0.02106463
## + Calcium           0.02096101
## <none>              0.00000000
## 
##                  Df    AIC      F Pr(>F)   
## + Organic_Matter  1 159.43 5.8777  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.08664823 
## Call: abund_distnifH ~ Organic_Matter 
##  
##                    R2.adjusted
## <All variables>     0.42412042
## + Barium            0.14270722
## + Av_Potassium      0.14238286
## + Av_Phosphorus     0.13935754
## + Nitrate_Nitrite   0.13484559
## + Zinc              0.13050213
## + Magnesium         0.12895812
## + Phosphorus        0.12745386
## + Manganese         0.12724465
## + Copper            0.12206419
## + Iron              0.12132345
## + Cobalt            0.12126082
## + Ammonia           0.11844446
## + Nickel            0.11807834
## + Sulfur            0.11669887
## + Moisture_Content  0.11655745
## + Calcium           0.11156606
## + pH                0.10979577
## <none>              0.08664823
## 
##          Df    AIC      F Pr(>F)   
## + Barium  1 157.08 4.2814  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.1427072 
## Call: abund_distnifH ~ Organic_Matter + Barium 
##  
##                    R2.adjusted
## <All variables>      0.4241204
## + Av_Potassium       0.2008529
## + Av_Phosphorus      0.1900137
## + Manganese          0.1880846
## + Phosphorus         0.1853480
## + Magnesium          0.1842902
## + Cobalt             0.1825912
## + Copper             0.1820141
## + Zinc               0.1809508
## + Iron               0.1802010
## + Ammonia            0.1754906
## + Moisture_Content   0.1750410
## + Nickel             0.1717581
## + Sulfur             0.1707810
## + Calcium            0.1678310
## + pH                 0.1675754
## + Nitrate_Nitrite    0.1669577
## <none>               0.1427072
## 
##                Df    AIC      F Pr(>F)   
## + Av_Potassium  1 154.35 4.5718  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.2008529 
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium 
##  
##                    R2.adjusted
## <All variables>      0.4241204
## + Manganese          0.2571847
## + Av_Phosphorus      0.2465866
## + Cobalt             0.2426324
## + Copper             0.2399983
## + Magnesium          0.2390988
## + Iron               0.2386871
## + Ammonia            0.2363363
## + Zinc               0.2360117
## + Nickel             0.2320877
## + pH                 0.2311367
## + Sulfur             0.2301444
## + Calcium            0.2287892
## + Phosphorus         0.2236127
## + Nitrate_Nitrite    0.2227016
## + Moisture_Content   0.2141867
## <none>               0.2008529
## 
##             Df    AIC     F Pr(>F)   
## + Manganese  1 151.47 4.637  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.2571847 
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese 
##  
##                    R2.adjusted
## <All variables>      0.4241204
## + Av_Phosphorus      0.3048552
## + Copper             0.2999316
## + Magnesium          0.2990454
## + Iron               0.2975378
## + Zinc               0.2954468
## + Cobalt             0.2944697
## + pH                 0.2893062
## + Nickel             0.2881294
## + Ammonia            0.2861456
## + Sulfur             0.2846872
## + Calcium            0.2815357
## + Phosphorus         0.2812081
## + Nitrate_Nitrite    0.2764049
## + Moisture_Content   0.2720061
## <none>               0.2571847
## 
##                 Df    AIC      F Pr(>F)   
## + Av_Phosphorus  1 148.93 4.2071  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.3048552 
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese +      Av_Phosphorus 
##  
##                    R2.adjusted
## <All variables>      0.4241204
## + Iron               0.3478017
## + Cobalt             0.3450564
## + Copper             0.3442042
## + Magnesium          0.3418955
## + Zinc               0.3394628
## + pH                 0.3331523
## + Sulfur             0.3315043
## + Phosphorus         0.3300620
## + Nickel             0.3283710
## + Calcium            0.3275764
## + Nitrate_Nitrite    0.3265021
## + Moisture_Content   0.3154696
## + Ammonia            0.3126530
## <none>               0.3048552
## 
##        Df    AIC      F Pr(>F)   
## + Iron  1 146.52 4.0021  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.3478017 
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese +      Av_Phosphorus + Iron 
##  
##                    R2.adjusted
## <All variables>      0.4241204
## + Sulfur             0.3885315
## + Calcium            0.3839583
## + Nickel             0.3693593
## + pH                 0.3689149
## + Nitrate_Nitrite    0.3651795
## + Zinc               0.3638021
## + Phosphorus         0.3632509
## + Moisture_Content   0.3622189
## + Magnesium          0.3594648
## + Ammonia            0.3562554
## + Copper             0.3550172
## + Cobalt             0.3480682
## <none>               0.3478017
## 
##          Df    AIC      F Pr(>F)   
## + Sulfur  1 144.06 3.9597  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.3885315 
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese +      Av_Phosphorus + Iron + Sulfur 
##  
##                    R2.adjusted
## <All variables>      0.4241204
## + Magnesium          0.4112957
## + Zinc               0.4085774
## + Nitrate_Nitrite    0.4069902
## + pH                 0.4051098
## + Nickel             0.4019619
## + Phosphorus         0.4015875
## + Copper             0.3991158
## + Moisture_Content   0.3920267
## + Ammonia            0.3894283
## <none>               0.3885315
## + Calcium            0.3869749
## + Cobalt             0.3831228
## 
##             Df    AIC      F Pr(>F)   
## + Magnesium  1 142.96 2.6569  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.4112957 
## Call: abund_distnifH ~ Organic_Matter + Barium + Av_Potassium + Manganese +      Av_Phosphorus + Iron + Sulfur + Magnesium 
##  
##                    R2.adjusted
## + Phosphorus         0.4318158
## + Nitrate_Nitrite    0.4268994
## <All variables>      0.4241204
## + Nickel             0.4189326
## + pH                 0.4174772
## + Zinc               0.4166306
## + Ammonia            0.4127613
## <none>               0.4112957
## + Moisture_Content   0.4108080
## + Calcium            0.4093634
## + Copper             0.4068650
## + Cobalt             0.4064199

Constrained Ordination

http://deneflab.github.io/MicrobeMiseq/demos/mothur_2_phyloseq.html#constrained_ordinations

CCA with selected model from Chemical CCA model

# CCA ordinate
cca_ord_nifH <- ordinate(
    physeq = physeq_nifH_ord, 
    method = "CCA",
    distance = abund_distnifH,
    formula = ~  Organic_Matter + Nitrate_Nitrite + Av_Potassium + 
      Av_Phosphorus + Manganese + Barium + Sulfur + Magnesium  )

# CCA plot
cca_plot_nifH <- plot_ordination(
  physeq = physeq_nifH_ord, 
  ordination = cca_ord_nifH, 
    color = "Site",
    axes = c(1,2)) + 
    geom_point(aes(colour = Site), size = 3) + 
    scale_color_manual(values =  farm_col_paired) 

# Now add the environmental variables as arrows
arrowmat_nifH_cca <- vegan::scores(cca_ord_nifH, display = "bp")

#Get appropiate scalling multipler 
mul<-vegan::ordiArrowMul(arrowmat_nifH_cca)

#Multiply biplot by scaling multiplier
arrowmat_nifH_cca_scale<-arrowmat_nifH_cca*3

# Add labels, make a data.frame
arrowdf_nifH_cca <- data.frame(labels = rownames(arrowmat_nifH_cca_scale), arrowmat_nifH_cca_scale)

# Define the arrow aesthetic mapping
arrow_map <- aes(xend = CCA1, 
    yend = CCA2, 
    x = 0, 
    y = 0, 
    shape = NULL, 
    color = NULL, 
    label = labels)

label_map <- aes(x = 1.3* CCA1, 
    y = 1.3 * CCA2, 
    shape = NULL, 
    color = NULL, 
    label = labels)

arrowhead = arrow(length = unit(0.02, "npc"))

# Make a new graphic
cca_plot_nifH + 
  geom_segment(
    mapping = arrow_map, 
    size = .5, 
    data = arrowdf_nifH_cca, 
    color = "gray", 
    arrow = arrowhead
  ) + 
  geom_text(
    mapping = label_map, 
    size = 4,  
    data = arrowdf_nifH_cca, 
    show.legend = FALSE
  )+
  ggtitle("CCA plot constrained ordination of nifH with selected Chemistry Model")+
  theme_bw()
## Warning: Ignoring unknown aesthetics: label

RDA with selected model from Chemical RDA model

# RDA ordinate
rda_ord_nifH <- ordinate(
    physeq = physeq_nifH_ord, 
    method = "RDA",
    distance = abund_distnifH,
    formula = ~   Av_Phosphorus + Sulfur + Magnesium + Zinc + Av_Potassium + Manganese )

# RDA plot
rda_plot_nifH <- plot_ordination(
  physeq = physeq_nifH_ord, 
  ordination = rda_ord_nifH, 
    color = "Site",
    axes = c(1,2)) + 
    geom_point(aes(colour = Site), size = 3) + 
    scale_color_manual(values =  farm_col_paired) 

# Now add the environmental variables as arrows
arrowmat_nifH_rda <- vegan::scores(rda_ord_nifH, display = "bp")

#Get appropiate scalling multipler 
mul<-vegan::ordiArrowMul(arrowmat_nifH_rda)

#Multiply biplot by scaling multiplier
arrowmat_nifH_rda_scale<-arrowmat_nifH_rda*700

# Add labels, make a data.frame
arrowdf_nifH_rda <- data.frame(labels = rownames(arrowmat_nifH_rda_scale), arrowmat_nifH_rda_scale)

# Define the arrow aesthetic mapping
arrow_map <- aes(xend = RDA1, 
    yend = RDA2, 
    x = 0, 
    y = 0, 
    shape = NULL, 
    color = NULL, 
    label = labels)

label_map <- aes(x = 1.3* RDA1, 
    y = 1.3 * RDA2, 
    shape = NULL, 
    color = NULL, 
    label = labels)

arrowhead = arrow(length = unit(0.02, "npc"))

# Make a new graphic
rda_plot_nifH + 
  geom_segment(
    mapping = arrow_map, 
    size = .5, 
    data = arrowdf_nifH_rda, 
    color = "gray", 
    arrow = arrowhead
  ) + 
  geom_text(
    mapping = label_map, 
    size = 4,  
    data = arrowdf_nifH_rda, 
    show.legend = FALSE
  )+
  ggtitle("RDA plot constrained ordination of nifH with selected Chemistry Model")+
  theme_bw()
## Warning: Ignoring unknown aesthetics: label

# CAP ordinate
cap_ord_nifH <- ordinate(
    physeq = physeq_nifH_ord, 
    method = "CAP",
    distance = abund_distnifH,
    formula = ~ Organic_Matter + Barium + Av_Potassium +
      Manganese + Av_Phosphorus + Iron + Sulfur + Magnesium  )

# CAP plot
cap_plot_nifH <- plot_ordination(
  physeq = physeq_nifH_ord, 
  ordination = cap_ord_nifH, 
    color = "Site", 
    axes = c(1,2)) + 
    geom_point(aes(colour = Site), size = 3) + 
    scale_color_manual(values =  farm_col_paired) 

# Now add the environmental variables as arrows
arrowmat_nifH_cap <- vegan::scores(cap_ord_nifH, display = "bp")

#Get appropiate scalling multipler 
mul<-vegan::ordiArrowMul(arrowmat_nifH_cap)

#Multiply biplot by scaling multiplier
arrowmat_nifH_cap_scale<-arrowmat_nifH_cap*1.9

# Add labels, make a data.frame
arrowdf_nifH_cap <- data.frame(labels = rownames(arrowmat_nifH_cap_scale), arrowmat_nifH_cap_scale)

# Define the arrow aesthetic mapping
arrow_map <- aes(xend = CAP1, 
    yend = CAP2, 
    x = 0, 
    y = 0, 
    shape = NULL, 
    color = NULL, 
    label = labels)

label_map <- aes(x = 1.3 * CAP1, 
    y = 1.3 * CAP2, 
    shape = NULL, 
    color = NULL, 
    label = labels)

arrowhead = arrow(length = unit(0.02, "npc"))

# Make a new graphic
cap_plot_nifH + 
  geom_segment(
    mapping = arrow_map, 
    size = .5, 
    data = arrowdf_nifH_cap, 
    color = "gray", 
    arrow = arrowhead
  ) + 
  geom_text(
    mapping = label_map, 
    size = 4,  
    data = arrowdf_nifH_cap, 
    show.legend = FALSE
  )+
  ggtitle("CAP Plot Constrained Ordination of nifH with Selected Model")+
  theme_bw()
## Warning: Ignoring unknown aesthetics: label